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We report the results of Molecular Dynamics (MD) simulations and formal modeling of the free 
energy surfaces and reaction rates of primary charge separation in the reaction center of Rhodobacter 
sphaeroides. Two simulation protocols were used to produce MD trajectories. Standard force 
field potentials were employed in the first protocol. In the second protocol, the special pair was 
made polarizable to reproduce a high polarizability of its photoexcited state observed by Stark 
spectroscopy. The charge distribution between covalent and charge-transfer states of the special 
pair was dynamically adjusted during the simulation run. We found from both protocols that the 
breadth of electrostatic fluctuations of the protein/water environment far exceeds previous estimates 
resulting in about 1.6 eV reorganization energy of electron transfer in the first protocol and 2.5 eV 
in the second protocol. Most of these electrostatic fiuctuations become dynamically frozen on 
the time-scale of primary charge separation resulting in much smaller solvation contributions to 
the activation barrier. While water dominates solvation thermodynamics on long observation times, 
protein emerges as the major thermal bath coupled to electron transfer on the picosecond time of the 
reaction. Marcus parabolas were obtained for the free energy surfaces of electron transfer by using 
the first protocol while a highly asymmetric surface was obtained in the second protocol. A non- 
ergodic formulation of the diffusion-reaction electron transfer kinetics has allowed us to reproduce 
the experimental results for both the temperature dependence of the rate and the non-exponential 
decay of the population of the photoexcited special pair. 
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I. INTRODUCTION 

The problem of bacterial photosynthesis has received 
enormous attention from both experimental and theoret- 
ical communitiesji'^i^'^'^'^ Here, we consider only the first 
step in the sequence of electronic transitions following the 
absorption of a visible photon by the special pair of the 
reaction center, the primary charge separation. The cal- 
culation of the rate of primary charge separation, which 
over several decades of intense research has effectively 
become the hydrogen molecule of bioenergetics, involves 
two components: the electronic communication between 
the primary donor and acceptor responsible for electron 
tunnelingii^i^ and the Franck-Condon factor describing 
the probability of bringing the donor and acceptor levels 
into resonance with each othcrJ^iii Our paper is con- 
cerned with that latter part of the problem which we dub 
as the energetics of primary charge separationj ^^i^'^i^"'' 

In addressing the issue of the energetics of charge sep- 
aration, we first want to dissect this complex problem 
into two, not necessarily simpler, questions: (1) What is 
the importance of the structural arrangement of the co- 
factors in the reaction center protein? and (2) What are 
the roles played by the protein and hydrating water in ac- 
tivating electronic transitions? Each of these questions 
has generated a significant amount of literature on its 
own, and we will not be able to provide a comprehensive 
discussion of each topic, focusing instead on our main 
goal, the factors affecting the free energy of activation. 

Since optical spectroscopy of bacteriochlorophyll co- 
factors can be studied separately, the most intriguing 



question related to our discussion is how the energetics 
of optical transitions and electron transfer are affected 
when the cofactors are assembled within the protein ma- 
trix. The notion often circulated in the literature^ is that 
protein provides a low-polarity environment lowering the 
free energy of solvation of embedded cofactors compared 
to solvation in water. Even though this statement is 
generally correct, we will show below that nuclear sol- 
vation approaching the thermodynamic limit of infinite 
observation (waiting) time is still quite significant for the 
electron transfer dipole formed by difference occupation 
numbers (atomic charges) of the electron in the donor 
and acceptor states. In particular, solvation of the elec- 
tron transfer dipole by water is not fully screened by the 
protein and still makes about 1 eV. In addition, the pro- 
tein matrix cannot be really considered non-polar since 
there is a significant contribution to the reorganization 
energy from the nuclear modes of the protein. It turns 
out that the notion of weak nuclear solvation of primary 
charge separation, required to explain the observed rates, 
cannot fully rest on the thermodynamic arguments, and 
the dynamics of the protein/ water thermal bath need to 
be involved. 

Solvation dynamics of optical chromophores in dense 
molecular solvents have been actively studied in the past 
decades ji^ii^iii The basic picture, first discovered in nu- 
merical simulations^^ and later confirmed by laboratory 
measurements^^ is that the decay of the solvation cor- 
relation function (Stokes shift correlation function, S{t)) 
involves two major components. The fast Gaussian com- 
ponent is caused by ballistic motions of the solvent in 
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FIG. 1: Relaxation times (a,c) and the Stokes shift dynamics 
(b,d) of the structural glass-formers (a,b) and proteins (c,d). 
The vertical lines marked "exp." denote the temperature at 
which the Stokes shift correlation function is recorded. The 
dashed lines in (b,d) show the fast Gaussian decay of S{t). 



the first solvation shell of the solute (quasi-loealized vi- 
brations in the case of a protein). The slow tail of S{t) 
is related to coUeetive a relaxation mostly caused by re- 
laxation of orientations of molecular permanent dipoles 
(dielectric relaxation) and quadrupoles.— The notion of 
a relaxation, that is the slowest relaxation on the micro- 
scopic scale, is not commonly invoked in the discussion of 
high-temperature solvation dynamics of small molecular 
dyesr^ but becomes critical in building a conceptual ba- 
sis for understanding the solvation dynamics of cofactors 
assembled within the hydrated proteinics 

Phenomenology developed for structural glass- 
formers^i helps to formulate the problem we are dealing 
with here. The typical temperature dependence of the 
relaxation time of a polar molecular liquid is shown in 
Figure [IJi. A high-temperature liquid has two relaxation 
times: reorientations of molecular permanent dipoles 
resulting in slow a relaxation and fast f3f relaxation 
related to collective anharmonic cage rattling. Corre- 
spondingly, the Stokes shift correlation function has 
two components: fast Gaussian decay coupled to f3f 
molecular motions and a slow tail coupled to a motions. 
This latter component is often connected to dielectric 
relaxation of the homogeneous solvent When the 
liquid is supercooled, the a component, which often 
becomes non-Arrhenius, separates from the slow /3 relax- 
ation {(3s) characterized by the Arrhenius temperature 
dependence^ (Fig. W^)- If all the components of the 
Stokes shift correlation function could be resolved at 
that low temperature, three major parts, corresponding 
to a, Ps, and Pf relaxation could have been seen, ft 
is this imaginary experiment, which is hard to realize 
in molecular liquids that bears a close connection to 
charge-transfer dynamics in proteins. 

For proteins, as well as for most polymer glass-formers, 
a and (3 relaxation are well separated in the temper- 
ature range of protein stabilityj^iSiiS In addition, the 
secondary P relaxation is typically split into several com- 
ponents with increasingly faster dynamics accompanied 
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FIG. 2: Schematic arrangement of cofactors in the bacte- 
rial reaction center. P is the special pair, B and H are 
monomeric bacteriochlorophylls and bacteriopheophytins, re- 
spectively. Electron transfer in wild-type reaction centers oc- 
curs almost exclusively along the L-branch of cofactors (sub- 
script "L"), while the M-branch (subscript "M") is mostly 
inactive. 



by smaller activation barriers (Figure [T]^). The rugged 
surface of the protein also complicates the dynamics, and 
a relaxation is known to disappear from the response of 
water in nano-confinementi^ The coupling of the trans- 
ferred electron to different modes of the protein/ water 
solvent may vary, and it is a priori not clear which mode 
will dominate the solvation dynamics. However, one can 
clearly expect Stokes shift dynamics to show at least 
three components including a Gaussian decay due to /?/ 
modes, some subset of Ps modes, and an a relaxation 
(Figure [TJi). The relative relaxation times and weights of 
these modes in the overall Stokes shift correlation func- 
tion are critical for the energetics of charge transfer as 
we show below. 

The geometric arrangement of cofactors in the mem- 
brane protein of the reaction center has been consid- 
ered in the literature mostly from the perspective of cal- 
culating the probability of electron tunneling incorpo- 
rated into the electron-transfer matrix element ii Early 
studies considered the possibility of direct charge sepa- 
ration from the special pair (P) to bacteriopheophytin 
(Hi) of the L branch of monomeric chromophores via 
a super-exchange mechanism involving nearby bacteri- 
ochlorophyll (Bl)'^^ More recent studiesi^'^^'^'^'^ 
have identified as an intermediate state in the se- 
quence of electron hops^^C a slower process from P* to 
Bl followed by a faster transition from Bl to Hl- The 
ener gy level of was placed between 331-450 cm~^ 
(refslIllO) and 650-800 cm"! cm"! (ref IH) below the 
energy level of the excited special pair P*, favoring in 
both cases sequential over superexchange transfer. In 
the present study, we will restrict our attention to the 
first of two hops limiting our calculations to the rate of 
transition from P* to Bl (Figure [2]). 
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The role of the structural arrangement of the special 
pair in the energetics of primary charge separation has at- 
tracted relatively little attention (see, however, Warshel's 
work^). The spectroscopy of the P^P* excitation and 
of the primary pair cation radical have been intensely 
studied)^i2^ along with extensive modeling of the energy 
transfer within the antenna complex and to the special 
pairi^ The question we address here is somewhat dif- 
ferent. Given that the special pair has evolved within 
the reaction center, we wonder if its particular structural 
arrangement makes any significant impact on the acti- 
vation barrier of primary charge separation. Since the 
sandwich of two bacteriochlorophylls making P is highly 
conserved in bacterial and plant photosynthesis^^ it might 
have some other role in the functionality of the reaction 
center aside from capturing the excitation from the an- 
tenna complex. 

The motivation for posing this question is provided 
by Stark experiments by Boxer and co-workers who 
showed a dramatic increase of the polarizability of P 
upon photoexcitationi^^ In fact, the polarizability change 
of about 10^ A3 upon photocxcitation^ places the spe- 
cial pair among the most polarizable molecules known 
(carotenoids, also present in the reaction center, make an- 
other group of champions). This remarkable observation 
is combined here with our previous studies of electron 
transfer in polarizable donor-acceptor complexes r^iiSiii 
which showed that the change in polarizability accom- 
panying charge transfer results in asymmetric, non- 
parabolic free energy surfaces for electron transfer. De- 
termining whether this polarization asymmetry can sig- 
nificantly effect the activation barrier is one of the goals 
of this study. 

In summary, by combining extensive Molecular Dy- 
namics (MD) simulations with formal modeling, we want 
to establish the basic ingredients contributing to the ac- 
tivation barrier of primary charge separation. The ques- 
tions we address are the following: (1) What is the set 
of primary nuclear modes (either protein or aqueous wa- 
ter/protein interface) that promote transfer of an elec- 
tron? (2) How to describe the activation events happen- 
ing on such a short reaction time? In particular, we show 
that non-ergodic chemical kinetics is required in this case 
to replace the standard Marcus picture based on equilib- 
rium distributions. (3) What is the effect of high polariz- 
ability of the photoinduced special pair on the energetics 
of the transition? 

We address the questions posed above by incorporat- 
ing the Stokes shift dynamics from MD simulations into a 
formal theory which we describe first below. The results 
of the calculations presented next are tested for consis- 
tency against experimental data. We use our simulation 
data obtained at different temperatures to compare the 
calculated rates with the results of Fleming et alX In 
addition, the recently published data by Wang et alM 
for the population decay of the photoexcited state of the 
special pair in a number of mutants offer an opportunity 
to use the dynamical electron transfer models^i^i^ to 




FIG. 3: Reaction complex of Rhodobacter sphaeroidest— The 
protein (gray) is surrounded by the micelle of LDAO deter- 
gent molecules (purple). The electron is transferred in se- 
quence from the photoexcited special pair (P, gold) to bac- 
teriochlorophyll (Bl, red) followed by even faster transfer to 
bacteriopheophytin (Hl, green). In MD simulations, the re- 
action complex is surrounded by 6 sodium ions, 30 NaCl pairs, 
and 10,506 water molecules which are not shown here. 

study the multicxponential population decay. These ex- 
perimental results arc also analyzed here by combining 
the Stokes shift dynamics from MD simulations with a 
formal model of non-ergodic chemical kinetics. The pic- 
ture that has emerged from all this effort is summarized 
in the discussion section of this paper. 

II. BASICS OF THE FORMALISM 

We approach the problem of calculating the rates of 
charge separation by combining the input from MD sim- 
ulations with analytical formalism. Simulations of the re- 
action center of Rhodobacter sphaeroides^ were carried 
out using Amber 8.04i We have followed the procedure 
first suggested by Ceccarelli and Marchi^ in which the 
reaction center is surrounded by the micelle of detergent 
(lauryl dimethyl amino oxide, LDAO) molecules mimick- 
ing the hydrophobic membrane, and also closely match- 
ing the experimental setup^ for photochemical studies of 
bacterial photosynthesis. The structure of the reaction 
center surrounded by LDAO molecules is shown in Fig- 
urcO The details of the simulation protocol are provided 
in Appendix |^ and the charging scheme of the cofactors 
and the protein/ water solvent is outlined in Appendix IB] 

A. Energetics of Primary Charge Separation 

Electron transfer is a tunneling event realized, in the 
Born-Oppenheimer approximation, at the resonance of 
the electronic donor and acceptor energies. The gap 
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between the electronic energies of the donor and accep- 
tor states, AE, makes the one-dimensional reaction co- 
ordinate X = AE that incorporates the whole mani- 
fold of possible nuclear modes affecting the electronic 
transition)^ Because many nuclear motions affect the 
donor and acceptor states in dense condensed media, the 
fluctuations of the stochastic variable AE are often well 
represented by the Gaussian statisticsi^ Therefore, the 
probability of reaching zero energy gap AE = 0, when 
electron tunneling takes place, is given by the Gaussian 
probability 

PiAE = 0) oc e-<^^>'/2C^-(0) (1) 

Here, the variance Cx{0) = ((<5X)2), SX = AE - (AE) 
is equal to the t = value of the time self-correlation 
function of the energy gap X{t) = AE{t) 

Cx{t) = {5X{t)5X{Q)) (2) 

and the brackets stand for an ensemble average. 

If one multiplies the probability of reaching the res- 
onance with the frequency ujf, of electron tunneling be- 
tween the donor and acceptor electronic levels, one ar- 
rives at the Marcus-Levich equation for the electron 
transfer rate^ 

fcET=C.ee-<^^>'/2^-(°) (3) 

In the case of non-adiabatic electron transfer considered 
here (weak electronic coupling between the donor and 
acceptor), the electronic tunneling frequency is given in 
terms of the electron transfer matrix element V by the 
following equation: 

- (4) 

Equation [3] is quite general and is limited only by the 
assumption of Gaussian fluctuations of the energy gap. 
In order to make it practical, one needs to connect the av- 
erage energy gap (AE') and the variance Cx(0) to physi- 
cal interactions present in the system made of the donor- 
acceptor complex and a thermal bath of nuclear degrees 
of freedom coupled to the transferred electron. Despite 
the obvious complexity of the system, a generally applica- 
ble separation of the average energy gap into three com- 
ponents is possible: the gas-phase gap AE^^^, the shift 
by non-polar interaction potentials AE"""^, and the shift 
by Coulomb interactions between the permanent partial 
charges of the solute and the solvent, AE'^ 

(AE) = AE^"^ + AE'""^ + AE^ (5) 

The gas-phase energy gap AE^^ is the difference be- 
tween the ionization potential of the donor and the elec- 
tron affinity of the acceptor in the gas phase. The 
two other components represent the interaction with the 
protein/ water solvent and can thus be combined into a 
solvent-induced (subscript "s") shift 



The separation of the average energy gap into a non- 
polar and Coulomb part is in fact related to the sepa- 
ration of time-scales first discussed in early work on po- 
larons in solids by Pekar^^^i^ Frohlich)^ and Feynmaui^ 
The fast electronic degrees of the solvent (which is com- 
posed of protein, detergent, and water in our problem) 
result in instantaneous equilibration of the transferred 
electron by induction and dispersion (London) forces. 
For our present application, the former is more signif- 
icant (superscript "ind" in eq [5] and throughout below) 
and we therefore explicitly consider this component. The 
last term in eqjSJ related to Coulomb interactions, fluc- 
tuates due to slow molecular motions of molecular rota- 
tions and translations. This term is often described in the 
electron transfer literature by the coupling of the electric 
field of the donor-acceptor complex to the inertial dipo- 
lar polarization^ which we consider after the induction 
component. 

The induction forces are produced by polarizing the 
medium by the electric field of the donor-acceptor com- 
plex. If atoms and/or molecules of the medium carry 
polarizabilities aj, the induction energy is the sum of 
polarization free energies of all such polarizable groups 
located at positions rj. The induction shift of the aver- 
age energy gap is then given by the change in the polar- 
ization free energy caused by changing the electric field 
of the donor-acceptor complex 

AE'"^ = - [^02(r,) - i?o\(r,)] ^ (7) 

This component of the energy gap is often not given ad- 
equate attention in the electron transfer literature, even 
though it can be quite significan t^^'^^ as we show be- 
low. The induction shift also depends on temperature for 
constant-pressure experiments because of thermal expan- 
sion, and this fact needs to be included in the modeling 
of temperature-dependent reaction rates. Even though 
the induction potential is established instantaneously by 
induced electronic dipoles, the interaction energy is mod- 
ulated by nuclear motions of the solvent (water and pro- 
tein) producing a non-zero component in the Gaussian 
distribution width (see below). 

The Coulomb part of the donor-acceptor energy gap 
has received the most attention over several decades of 
the theory development, and we will briefly set up the 
stage for our treatment of this part of the problem here. 
The linear response approximation, either in terms of 
the electrostatic interaction with the medium dipolar 
polarization^ or in terms of partial atomic charges,— 
has mostly been used as the basis for theory develop- 
ment. In the former case, one considers the polarization 
of the solvent by the electric field of the solute Egi in the 
initial electron transfer state. This equilibrium polariza- 
tion Pcq(r) at point r within the solvent is connected 
to the field Eoi(r') at point r' by generally a nonlocal 
response function x(r, r^) 60, 61,62 



AEs = AE""^ + AE' 



(6) 



Peq(r) =x(r,r')*Eoi(r') 



(8) 
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where the asterisk denotes volume integration over the 
variable r' and tensor contraction over the Cartesian 
components of the field with the corresponding compo- 
nents of the 2-rank tensor x- 

Equation [5] for the solvent polarization induced by the 
solute typically appears in theories of linear solvation in 
homogeneous liquids. In contrast, the protein matrix it- 
self and the protein/water interface are inhomogcneous 
with a possibility of generating a polarization field P^q 
unrelated to the electric field of the cofactors. This polar- 
ization will create an additional inhomogcneous compo- 
nent of the vertical shift Ai^j'^j^ that cannot be calculated 
from the linear response approximation. The incrtial (nu- 
clear) polarization field Pg^ -|- Poq docs not change on the 
time-scale of electronic transition and creates a shift of 
the donor-acceptor energy gap by the amount determined 
by the change in electron's electric field AEq upon the 
transition: 

AE^ = - (P^q + Peq) * AEo - AEgi^ - Eoi * X * AEq 

(9) 

In the original Marcus formulation^ the average ver- 
tical energy gap was separated into the Coulomb reorga- 
nization energy and the Coulomb part of the Gibbs 
energy of the reaction, AG'~'. By using the identity 
Eoi = Eo - AEo/2, Eo = (Eqi + Eo2)/2 in cq [i one 
gets 

AE'^ = + AG^ (10) 

where 

A'^ = (l/2)AEo*x* AEo (H) 

and 

AG^ - AEg^ - Eo * X * AEo (12) 

The Coulomb part of the Gibbs energy then combines 
with the gas-phase gap and the induction shift to make 
the overall reaction Gibbs energy 

AG = A£;s'^^ + AS'"'^ + AG^ (13) 

Combined together, eqs[5l[Tni and [Ql lead to the stan- 
dard Marcus relation for the vertical average energy gap 

{AE)^AG + \^ (14) 

The separation of the average energy gap into the equi- 
librium Gibbs energy and reorganization energy compo- 
nents makes sense when the former can be measured sep- 
arately. In spectroscopy, the average gap is given by the 
maximum of the corresponding spectroscopic band (or, 
more precisely, by the first spectral moment) and this 
separation is not necessary. Likewise, the average energy 
gap is directly accessible from MD simulations, so the for- 
mulation in terms of the average energy gap is also more 
convenient from the simulation perspective. Even more 
importantly, the Gibbs energy of the reaction loses its 



direct connection to equilibrium thermodynamics in non- 
ergodic reaction kinetics, which we formulate and apply 
to the calculation of the rates below. In this framework, 
the formulation of electron transfer thermodynamics in 
terms of the first and second cumulants of the donor- 
acceptor energy gap is the only formal approach to the 
problem available at the moment. 

The Gaussian width, Cx{0) in eql^l generally needs a 
separate determination. It is calculated as the variance of 
the sum of all solute-solvent interaction potentials affect- 
ing the energy of the transferred electron. The problem is 
simplified for the Coulomb interactions. These are long- 
rangcd and arc typically well described by the linear re- 
sponse approximation. Therefore, the high-temperature 
limit of the fluctuation dissipation theorem^ applies to 
the Coulomb part on nuclear fluctuations with the result- 
ing factorization of Cx (0) into temperature and reorga- 
nization energy^2i^ 

Cx{0) = 2kBTXs (15) 

A significant simplification of this route is achieved 
through the fact that the variance is determined in terms 
of the same response function as the one used for the 
Coulomb part of the average energy gap (eq[5]), thus re- 
ducing the number of independent response functions re- 
quired by the theory. 

This procedure does not apply to short-range induction 
forces which do not follow the macroscopic fluctuation- 
dissipation theorem; the calculation of their first and sec- 
ond cumulants requires microscopic response functionsi^ 
The main consequence is that the induction compo- 
nent does not factorize into temperature and a weakly 
temperature-dependent energy parameter. The result is 
a generally non-Arrhenius form of the rate constant^ in 
eq|3]in which the variance can be written as 

Gjf(0)~2fcBTAC + G'"'^(0) (16) 

Since the Coulomb and induction interaction sum up in 
the energy gap, a cross term needs to be taken into ac- 
count, and we have included it into G'"'^(0) as follows 

G'"^(0) = {{dE''"^f) + 2{5E^5E''"^) (17) 

Despite these complications which take away the solid 
foundation behind factoring the variance into the temper- 
ature and energy components^^i^ we will follow the es- 
tablished tradition and define the solvent reorganization 
energy as as the sum of induction and Coulomb terms 
(c/. to eqEl) 

A, = A""' + A^ (18) 

where 

A'"'^ = G'"'^(0)/(2/fcBr) (19) 

The probability of electron transfer can be affected 
by intramolecular vibrations of the soluteJ^ These can 
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be added to the formahsm outhncd here by summing 
up probabihties of transitions between separate vibronic 
channels. These transitions are known to significantly 
affect the transition probability in the inverted region, 
(AE') < 0, but can be neglected for transitions in the 
normal region, (Ai?) > 0, considered here. An extension 
to the former case is well developed in the literatureii 
and does not pose fundamental difficulties. 



B. Stokes shift dynamics 

The characteristic timescales of nuclear fluctuations af- 
fecting charge transfer can be extracted from the time 
correlation function in eq[2l or from its normalized value 
known as the Stokes shift correlation function 



S{t)^Cx{t)/Cx{0) 



(20) 



As mentioned above, the typical shape of S{t) in com- 
plex condensed media includes a fast Gaussian compo- 
nent and a multi-exponential (or stretched-exponential) 
tail. A two-exponential tail is used to fit our simulation 
results with Cx (t) in the form 



Cx(i) =C'"'i(t) + CC(t) 



(21) 



where 



C^{t) = 2fcBT 



A 



(22) 

Here, tq is the relaxation time of the Gaussian decay 
and Ti and T2 are two exponential relaxation times. In 
addition, Ag and Ai are the corresponding reorganization 
energy components such that A^ = Ag -t- Af 



A?. 



C. Non-ergodic activation kinetics 



The arguments presented in sec III Al arc based on 
equilibrium statistical mechanics representing the com- 
ponents of the activation barrier as equilibrium (free) 
energies. This formulation in fact assumes a certain 
separation of time-scales, that is the time of the reac- 
tion Tet = rnust be much longer than all relaxation 
times (tg, n, etc.) of the nuclear modes coupled to the 
transferred electron. This assumption certainly breaks 
down for our problem combining the extremely short 
time of natural primary charge separation (ca. 3 ps) with 
the disperse relaxation spectrum of the protein/ water 
solvent li^i^Si^ What we face here is the obvious case of 
ergodicity breaking^ of the nuclear fluctuations involved 
in the reaction activation, which raises the question of 
how to approach the calculation of the reaction rates. 

The Stokes shift correlation function provides a con- 
sistent approach to formulate the kinetics of non-ergodic 
electron transfer. We first note that the equilibrium lin- 
ear response function introduced in Sec. Ill Al in the 
direct space domain, can be extended to the time domain 



to cover the time correlation functions of the energy gap 
fluctuations. The equilibrium ensemble average produc- 
ing the solvent response component of the average energy 
gap (c/. to eqin]) can then be given as a frequency integral 
of the Fourier transform x(w) 



Jo 



(23) 



This representation offers a systematic approach to cal- 
culating the non-ergodic solvent response. The integral 
in eq is over all possible frequencies of nuclear mo- 
tions, implying that all of them contribute to the aver- 
age. In fact, the time-scale of the reaction tet limits 
the frequency spectrum only by those frequencies that 
are higher than the rate of the reaction A;et- The non- 
ergodic energy gap thus becomes 



dojEoi * x(w) * AEq 



(24) 



Along the same lines, the non-ergodic reorganization en- 
ergy can be defined by using the same step- wise frequency 
flltcr: 

/•OO 

A^(fcET) / dwAEo * x(a') * AEo (25) 

An alternative representation is through the Fourier 
transform of the Stokes shift correlation function 



C^ic) = e-*C^(i)di/(27r) 



as follows^ 



A^(fcET)=/3/ C'^iu;)duj 



(26) 



(27) 



where (3 = l/{k^T). Equations [MH?7l suggest that 
A^(fcET) can be obtained from the Stokes shift correlation 
function calculated from MD trajectories while a formal 
theory is required for x(u;) to determine A£'^(fcET)- 

The notion that the parameters entering the activation 
barrier become functions of the electron transfer rate cre- 
ates the necessity to consider the calculation of the rate 
constant as a self-consistent problem given as the solution 
of the following equation: 

KET = (fcET)exp [-(Ai;(fcET))V2Cx(0,A:ET)] (28) 

Here, the rate-dependent energy gap can be re-written 
based on eq [M] as 

(Aii;(fcET)) - AE^-^ + AE'^^ + AEg,^ + fZ{kET)AEg 

(29) 

where, based on our simulations discussed below, we as- 
sume that the induction component of the shift does not 
involve slow relaxation and only Coulomb solvation gets 
cut off by breaking ergodicity. Accordingly, the Gaussian 
width in eq [55] takes the form 



CxiO, kET) = C'"'i(0) + 2fcBTAC(fcET) 



(30) 
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where 

A^(fcET) = /„'c(feT)AC (31) 

In eqs [221 and [211 we have introduced the parameters of 
non-ergodicity of nuclear fluctuations contributing to the 
vertical energy gap, and to the reorganization energy, 
f^^. The parameter f^^ can be readily calculated from 
eqsinilll and [23 

/4 = Ag/A^ + (2/7r) iX,/\^) cot-i (fcETT.) (32) 

i=l,2 

The procedure outlined above can be used to construct 
the free energy surfaces of electron transfer. The widely 
accepted definition of the free energy surfaces for elec- 
tron transfer follows the general procedure of defining the 
Landau functional^ in which the hypcrsurfacc X — AE 
generates the incomplete partition function 

e-P''^''^ j 5{AE-X)e-^"dT (33) 

In this expression, Ai? depends on all nuclear modes Qm 
n = 1, . . . , M in the system. In addition, H is the system 
Hamiltonian in the initial state of the electron-transfer 
system and dV is the element of phase space. 

In applications to processes happening on short time 
scales, one needs to generalize eq[331to exclude a subset 
of frequencies not contributing to the process: 

n.cj<A;ET 

. . 

In this equation, the product of delta functions eliminates 
the low-frequency modes from the partition function. 

D. Polarizability of the special pair 

There is a significant body of experimental^i^i^Sii^ 
and computationa l'^^1^^'^'^1^^'^^1^^ evidence of a strong 
mixing of covalent, (PL-Pj\f )*, and charge-transfer, (P^/- 
P^)*, states within the photoexcited special pair, where 
Pa/ and Pl are the M and L subunits of the special 
pair (Figures [Hand [3]). Although the average amount of 
charge transfer between two subunits is small)^ about 
0.1 of the electronic charge in the gas phase and 0.2 
in the reaction center, the fluctuations of the elec- 
trostatic potential of the protein/ water solvent create 
significant fluctuations of the extent of charge trans- 
fer. Correspondingly, the fluctuating population of 
the charge-transfer state (Pj^-P^)* creates fluctuating 
charges /S.Zj = nc^AZj at the atomic sites of the special 
pair {j runs over the atoms of P). In this representation, 
AZj are the difference of atomic charges of the special 
pair between the ionized excited state (P^-PJ^)* ^^^d the 
covalent state (Pl-Pm)*, and ncT is the population of 
the charge-transfer state. Physically, this redistribution 



of charge in response to an external electrostatic field im- 
plies that the special pair is polarizable with the instanta- 
neous induced dipole moment equal to pcT = '^ctA/^ct- 
Here, A/^ct is the dipole moment between the ionized 
and neutral states of the special pair. The importance of 
the induced dipole moment for electron transfer is that 
the instantaneous electron transfer dipole becomes modi- 
fied from the dipole /xet created by the set of permanent 
charges Aq^ {k runs over all atoms of the cofactors in- 
volved in electron transfer) to a new fluctuating dipole 
moment /xet + Pct- Since the solvent reorganization 
energy is proportional to the average squared dipole mo- 
ment 

(X (^{fiET + Pcrf) (35) 

the appearance of the induced dipole can potentially 
modify the energetics of electron transfer.— i^i^ 

In order to model the effect of polarizability of the 
special pair on the statistics of the donor-acceptor energy 
gap, we have adopted the following simulation algorithm. 
The charges Zj of the primary pair are re-calculated at 
each fifth MD step according to the equation 

Zj = z^ + uctAZj (36) 

where are the charges of two decoupled bacteri- 
ochlorophylls obtained from our DFT calculations (Ap- 
pendix [B1 and supporting information). The extent of 
charge dclocalization ncT is calculated by diagonaliz- 
ing, at each fifth step of the MD trajectory, the two- 
state quantum Hamiltonian characterized by the elec- 
tronic coupling J and the instantaneous energy gap be- 
tween two states 

Ae = Aes^" + Ae''^'^ + (1/2) ^ AZ^cjj^ (37) 

j 

Here, Ae^^^ is the gas-phase energy separation between 
the neutral and ionized states of P and (f)j is the electro- 
static potential of the surrounding protein/water solvent 
at the position of atomic charge j. Ae'"'* in eq [371 is 
the induction shift of the energy gap and the parameters 
Ae^'^^ and J are tabulated in Appendix [51 

E. Polarizable special pair and free energy surfaces 
of electron transfer 

The description of Coulomb solvation presented in 
sec III Al change significantly when the special pair is 
polarizablc^MO The main modification here is that the 
atomic charges and hence the electric field of the co- 
factors become a function of the solvent polarization P 
through the extent of charge dclocalization tt-ct- The 
electric field Eqi changes from the value commonly cal- 
culated from the vacuum charge distribution to Eoi[P]: 

Eoi^Eoi[P] (38) 



8 



A general solution for the free energies of electron trans- 
fer in this case has not been found so far, although an 
analytical theory can be formulated in the case of dipole 
solvation!^ Alternatively, the field Eqi [P] can be linearly 
expanded in the solvent polarization P around its equi- 
librium value 

Eoi[P] -Eo[Peq + Py + F.(5P (39) 

where F is a 2-rank tensor and SP = P — Poq — P^q. 

When the form of the field given by eq[39]is substituted 
into the standard Hamiltonian^ of the solute linearly 
coupled to the Gaussian field P, one gets 

H = -Eoi [Peq + P^q] * P + (1/2)5P * x^od * (40) 

where = ~ 2F is the new, modified linear re- 

sponsc function of the Gaussian polarization field renor- 
malized by the solute polarizability. Since the polariz- 
ability tensor F is generally different in the initial and 
final electronic states, the donor-acceptor energy gap be- 
comes a bilinear function of the Gaussian field P in 
contrast to the linear function used to derive the Mar- 
cus parabolas. The main consequence of that change 
is that the statistics of energy gap fluctuations become 
non-Gaussian. The free energy surface loses its parabolic 
shape predicted by Marcus theory and can instead be 
represented by the analytical results of the Q-modeli^ 

G{X) = a (^^J\{AE) -a\c -X\ - (41) 

Here, a > is the non-parabolicity parameter describing 
the deviation of the free energy surface from the parabolic 
shape. The limit a ^ oo recovers the Marcus barrier 
thermodynamics. 

III. RESULTS 

We believe that this paper reports the most extensive 
MD simulations on the bacterial reaction center following 
previous simulation efforts in this field i^'^'^^i^^'''^i^°'^-^i^^ 
The overall length of 118 ns of MD trajectories, of which 
100 ns were used for the production analysis, required 
39.8 CPU years. All simulations were done in parallel 
using 128 CPUs of ASU's HFC facility. The analysis of 
the simulations was performed by a parallel code devel- 
oped for this project that reads directly binary AMBER 
files. The analysis was run in parallel on 10 Opteron 
CPUs and required overall 4.8 years of CPU time. 

Two sets of simulations have been done. The first set, 
which we will label SI, was performed at six different 
temperatures. It employed the standard protocol of MD 
force fields with fixed atomic charges. The equilibrium 
MD trajectories were used to calculate the statistics of 
the donor-acceptor energy gap and the Stokes shift cor- 
relation functions. In this calculation, in addition to 
Coulomb interactions, induction solute-solvent interac- 
tions were computed. The atomic polarizabilities were 



taken from a modified Thole parametrizationi^ The in- 
duction potential was not a part of the simulation algo- 
rithm, thus assuming that the exploration of the phase 
space of the nuclear motions can be accomplished with 
the standard force fields. Since these force fields effec- 
tively incorporate polarizability in terms of permanent 
charges, in order to avoid double counting, the charges 
of the solvent (protein and water) were multiplied by 0.89 
in analyzing the data, following the convention adopted 
in the literature.— 

Six trajectories of SI protocol were produced for the 
initial state of the reaction complex, (P-B^)*, at differ- 
ent temperatures. The atomic partial charges calculated 
by us at the DFT level (Appendix |B|) were supplemented 
by the force-field parameters of bacteriochlorophyll de- 
veloped by Marchi and co-workers The atomic charges 
of the ground-state bacteriochlorophyll were used for the 
exited state of the primary pair assuming that photoex- 
citation does not greatly alter the charge distributioni^ 
One simulation trajectory at 300 K was produced for 
the charge-separated state P+-B^ corresponding to the 
first hop of the electron in the sequential mechanism. For 
this simulation, the positive charge of F"^ was distributed 
among the two cofactors of the special pair as described 
in Appendix [B] and the charge distribution of the bac- 
teriochlorophyll anion was calculated at the DFT level 
(supporting information). 

The second set of simulations, labeled as S2, required 
changing the standard MD protocol (see Appendix |B]). In 
these simulations, quantum polarizability of the special 
pair was accounted for by diagonalizing the 2x2 Hamil- 
tonian matrix of the charge-transfer state between the 
two parts of F at each fifth step of the MD trajectory, 
a procedure known in the literature as the empirical va- 
lence bond approach i^^iS^ The Hamiltonian diagonaliza- 
tion allows one to calculate the extent of charge transfer 
between two bacteriochlorophylls in F and dynamically 
adjust charges of the special pair. This simulation pro- 
tocol thus incorporates an extremely high polarizability 
of F* revealed by Stark spectroscopy measurements of 
Boxer and co-workers i^^^^ 



A. Energetics 

Two energy parameters are of main importance within 
the Gaussian picture of electron transfer activation (Mar- 
cus model). These are the average donor- acceptor energy 
gap and the energy gap variance (eq[3]). These parame- 
ters, obtained from MD simulations at different temper- 
atures, are listed in Tables U and [TTl The complete set 
of first cumulants from both SI and S2 simulations is re- 
ported in Table U The S2 entry in Table M is limited 
to 300 K since the second cumulants at other tempera- 
tures did not converge on the time-scale of the simulation 
trajectories. 

Since we are dealing here with a heterogeneous solvent 
composed of a protein matrix and aqueous environment. 
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TABLE I: Components of the average energy gap for primary charge separation (all energies are in eV). 



Protocol T/K A£;{;^' 












AEs 


SI 


77 


0.170 


-0.500 


-0.330 


-0.058 


-1.234 


-1. 


.292 


-1.623 




200 


0.144 


-0.728 


-0.584 


-0.057 


-1.199 


-1, 


.256 


-1.840 




250 


0.382 


-0.749 


-0.367 


-0.054 


-1.183 


-1, 


.237 


-1.604 




300" 


0.205 


-0.678 


-0.473 


-0.055 


-1.164 


-1, 


.219 


-1.691 




300' 


-0.365 


-1.278 


-1.643 


-0.082 


-1.071 


-1, 


.153 


-2.796 




350 


0.440 


-0.856 


-0.416 


-0.050 


-1.094 


-1, 


.144 


-1.559 




400 


0.093 


-0.330 


-0.237 


-0.075 


-0.853 


-0, 


.928 


-1.164 


S2 


250 


0.399 


-0.502 


-0.103 


-0.048 


-1.042 


-1, 


.090 


-1.193 




275 


0.404 


-0.535 


-0.131 


-0.051 


-1.065 


-1, 


.116 


-1.247 




300" 


0.310 


-0.632 


-0.323 


-0.052 


-1.076 


-1, 


.128 


-1.451 




325^= 


0.310 


-0.168 


0.141 


-0.020 


-0.836 


-0, 


.857 


-0.716 




350^^ 


0.347 


-0.594 


-0.247 


-0.053 


-0.987 


-1, 


.041 


-1.287 



"Obtained from 10 ns long MD trajectories, the unmarked data 
refer to 5 ns of simulations. 
''Data for the final charge transfer state P+— B^, 5 ns trajectory. 
'^Obtained from 15 ns long MD trajectories. 



TABLE II: Reorganization energies calculated from fluctuations of the energy gap feg I15p . All energies are in eV. 

Protocol T/K ^w"^ A^jrot A^^^*^ ^prot Aprot 



51 77 0.019 0.065 0.070 0.187 0.168 0.351 0.191 0.245 0.421 
200 0.001 0.057 0.062 0.756 0.182 0.845 0.756 0.251 0.903 
250 0.016 0.076 0.081 1.634 0.341 1.938 1.639 0.419 1.955 
300' 0.047 0.112 0.119 1.136 0.375 1.564 1.124 0.466 1.598 
300" 0.047 0.146 0.149 1.393 0.441 1.542 1.379 0.593 1.692 
350 0.110 0.187 0.191 0.948 0.644 1.508 0.944 0.853 1.508 
400 0.139 0.249 0.275 0.767 0.567 1.010 0.866 0.797 1.335 

52 300 0.481 0.682 0.697 1.439 0.735 1.839 1.454 1.385 2.513 



n^ind _j_ deviates slightly from As because of numerical uncer- 
tainties of averaging. 

'Obtained from 10 ns long MD trajectories, the unmarked data 
refer to 5 ns of simulations. 

'^S ns data for the final charge transfer state P^-B^ . 



the separation of these two first cumulants of the energy 
gap into the water and protein contributions provides 
mechanistic insights into the factors influencing electron 
transfer activation. In addition, we split the relevant en- 
ergies into contributions from non-polar and Coulomb 
interactions. Finally, the introduction of polarizability 
(charge fluctuations) of the special pair shifts relative 
weights of each component in the activation barrier and, 
more importantly, results in significant deviations from 
the Gaussian picture of Marcus parabolas. 

Figure |4] reports the distribution of Coulomb and in- 
duction components of the energy gap from simulations 
of both the non-polarizable and polarizable special pair. 
The Coulomb interactions have Gaussian statistics where 
the width is consistent with the reorganization energies 
listed in Table |TT1 The average shifts arising from water 
and the protein have opposite signs. Therefore, the po- 
larization of water by the protein matrix contributes to 
the destabilizing of the charge-transfer state, as was also 
observed by Parson et al^ On the contrary, the protein 
matrix makes the dominant contribution into stabilizing 
the charge-separated state. The induction shift of the 
average energy gap, arising primarily from the protein 



matrix (Ind(prot) in Figure|4|), is about twice larger than 
the Coulomb shift which largely cancels out between its 
protein and water contributions (Table |l| . On the con- 
trary, the width of the distribution of induction energies 
is small relative to the Coulomb interactions for non- 
polarizable (SI) simulations (in accord with assessment 
of analytical theories^), but grows significantly for the 
polarizable (S2) simulation protocol (Table . 

The splitting of the total self-correlation function 
Cx{0) into the individual protein (subscript "prot") and 
water (subscript "w" ) components requires an estimate 
of the cross-correlation term Aw.prot between the water 
and protein interaction potentials: 

As Aprot ^" Aw ^" Api'ot.w 

(42) 

This latter part turns out to be significantly smaller than 
the individual protein and water components, as can be 
inferred from the last three columns in Table Ull bv com- 
paring the total solvent reorganization energy As with 
the sum of the two components, Aprot + Aw. 

Notwithstanding such little attention paid in the 
electron-transfer literature to non-polar interactions, the 
induction shift is the main part of the solvent effect on 
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FIG. 4: Normalized distributions of components of the donor- 
acceptor energy gap in non-polarizable (SI, solid lines) and 
polarizable (S2, dashed lines) simulation protocols (T = 300 
K). Marked in the plot are the Coulomb interaction due to 
the protein (C(prot)) and water (C(w)) and the induction 
interaction (Ind(prot) for the protein and Ind(w) for water). 
SI and S2 mark the distributions of the total energy gap for 
the non-polarizable (SI) and polarizable (S2) special pair. 



FIG. 5: Components of the solvent reorganization energy 
from water and protein from MD simulations vs the obser- 
vation time Tobs. Points refer to MD data (S2 protocol, 300 
K) and the dotted lines indicate the fits to eq 1441 The in- 
set shows the initial portion of the plot; the relaxation times 
used to fit the MD results to eq|44]are: tr — 218 ps (total), 
tr = 390 ps (protein), and tr = 764 ps (water). 



the average energy gap of charge separation. Its value 
can be estimated from some simple arguments. If one 
assumes that atomic polarizabilitics are distributed with 
a constant density around the donor and acceptor, one 
arrives at a simple expression^ 

Here, Rd and Ra are the radii of the donor and acceptor 
and Rda is the distance between them. In addition, nprot 
is the refractive index of the protein matrix and e is the 
elementary charge. For the average refractive index of 
the reaction center— ripj-ot — 1.473 and the radius of 
the bacteriochlorophyll unit Ru = Ra = 5.6 A obtained 
from its vdW volume one gets AE'^^'^^ = -1.09 eV at the 
crystallographic distance Rda = 11-3 A. This number 
compares favorably with the induction shift of AE'pJ.^t ~ 
-1.16 cV from MD simulations at T = 300 K (Table [B 
SI protocol). 

The positive slope of the induction shift of the average 
energy gap is caused by the temperature expansion of 
the protein. Based on the data shown in Table |T] for SI 
simulation protocol, the logarithmic derivative of the in- 
duction shift with temperature, din AiJ^'^/dT, is within 
the limits (4 - 8) x 10""' K"^ According to eg H51 this 
derivative should be equal to thermal expansivity of the 
protein (Clausius-Mossotti equation). Indeed, the loga- 
rithmic slope of the induction shift agrees reasonably well 
with the reported^ expansion coefficients of proteins of 
the order of 8 x lO"'' K-\ 

Several of the MD simulation results reported here 
turned out to be quite surprising. Among the unex- 
pected findings arc quite large values of the solvent (pro- 
tein and water) reorganization energies, contrasting the 
commonly low values (ca. 0.1-0.2 eV) circulating in the 
literature. In particular, water is far from be- 
ing screened by the protein^S making the main portion 



of the energy gap variance in the SI protocol, and be- 
ing surpassed by the protein in the S2 protocol. In fact, 
the values of water reorganization energy found here are 
more typical of small redox couples in aqueous solution^ 
than of the often anticipated hydrophobic screening by 
the protein matrix. 

What is different in our simulations compared to previ- 
ously reported simulation data^i^>^ is the length of the 
simulation trajectory which has allowed us to push the 
numbers for the reorganization energies closer to their 
thermodynamic limit. Indeed, as is seen in Figurc[5l the 
reorganization energy as a function of the length of the 
simulation trajectory (observation time Tobs) levels off 
by the end of the 5-10 ns production run. However, this 
long-time reorganization energy is not relevant for the 
short-time charge-separation dynamics since a significant 
subset of nuclear modes gets dynamically arrested on the 
picosecond time-scale at which the reorganization energy 
as a function of Tobs starts to sharply decline (Figure [5]). 
The dependence of the total reorganization energy and 
its components on the observation window can be fitted 
to a one component Debye equation (cf. to eg I32p 

A(robs) (X cot "^(TflV Tobs) (44) 

with the effective relaxation time t/j. The fits shown by 
dotted lines in Figure [5] indicate that the system starts 
to lose ergodicity on the time-scale of several hundred 
picoseconds. 

The data in Figure [5] have been generated according 
to the following procedure. First, a trajectory of indi- 
vidual protein/solvent vertical energies is created from 
the sum of their respective Coulomb and induction com- 
ponents. Second, a smaller trajectory window of length 
Tobs is cut from the full MD trajectory. Third, the en- 
ergy gap variance is calculated on this smaller observa- 
tion window which is then moved along the entire tra- 
jectory. Each time the window is shifted, the variance is 
calculated with the average energy gap set to its average 
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FIG. 7: Induction (squares) and total protein (diamonds) re- 
organization energies from the present MD simulations and 
experimental mean square displacements of hydrogens of bac- 
teriorhodopsin scaled with the inverse temperature, {{Sx)'^) /T 
(circles). The experimental data were obtained by neutron 
scattering.— All parameters have been normalized to their 
corresponding values at 200 K. 



FIG. 6: Protein (a) and solvent (b) reorganization energies 
from MD simulations as functions of the observation time 
Tobs and temperature. The closed square indicates the result 
of ref obtained on a 40 ps observation window and the 
closed circle refers to the 20 ps window used in ref Isol. 



on that particular window. The individual variances are 
then averaged among the results from each sliding win- 
dow, and the average reorganization energy is reported 
as the A (robs) in Figure [Sj 

Reorganization energies calculated from this algorithm 
are plotted vs temperature in Figure [51 As in Figure [5l 
shortening the observation window lowers the reorgani- 
zation energy. On the 4 ps observation window, most of 
the multi-exponential Stokes shift relaxation is dynam- 
ically arrested (sec below) and only ballistic Gaussian 
relaxation from the Coulomb interactions and the mod- 
ulation of induction interactions by density fluctuations 
contribute to the reorganization energy. For this short 
observation time, the reorganization energy falls in the 
range of values commonly reported from fitting the ex- 
perimental reaction rateS] ^°i^^ In particular, our results 
from 50 ps observation window are consistent with the 
previous report by Parson et al^ using 40 ps of the sim- 
ulation trajectory (for Rp. viridis), while the result from 
20 ps simulations from Treutlein et ali^ is slightly below 
that value (closed points in FigureOj). We do not expect 
a close agreement here since our algorithm of sliding win- 
dow generally gives rise to higher reorganization energies 
than a single observation. 

The reorganization energy from the protein is an in- 
creasing function of temperature for all observation win- 
dows (Figure [5^). On the contrary, for water reorga- 
nization, the negative temperature slope expected from 
equilibrium statistical mechanics^^ is reverted by non- 
ergodicity to a positive one (Figure [BJd). This effect is 
caused by a temperature-depending unfreezing of the nu- 
clear modes when relaxation becomes faster with increas- 



ing temperature. The downward turnover of A(T) for 
the 1 ns observation window (upper curve in Figure [6)3) 
marks the return of the system to equilibrium statistics 
with the negative slope of A(T) also seen in our previous 
simulations of a small solute in SPC/E water 1^ 

The opposite temperature dependence of the protein 
and water reorganization energies at long observation 
windows points to a distinctly different character of the 
corresponding nuclear modes. While water molecules 
alter the donor-acceptor energy gap mostly by libra- 
tional/rotational motions typical of polar liquids, the 
protein nuclear modes are predominantly vibrational. 
The temperature dependence of Aprot seems to correlate 
well with the temperature dependence of atomic displace- 
ments of the protein matrix as is illustrated in Figure [7| 
where we show the better converged induction reorgani- 
zation energy along with the total protein reorganization 
energy. The temperature change of these reorganization 
energies is compared with mean square displacements of 
hydrogens in bacteriorhodopsin obtained from inelastic 
neutron scatteringi^i 

It is by now well established that protein vibrations 
start to deviate from the straight line of harmonic mo- 
tions at the transition temperature of about ~ 200 — 
220 K marking the rise of anharmonic motions (including 
side-chain rotations).— i^^i^^ Therefore, the mean-square 
displacement scaled with inverse temperature, {{Sx)'^)/T, 
is a flat function at low temperatures starting to rise 
above the transition temperature T^. The same trend is 
seen for the total protein reorganization energy and its in- 
duction component, which both turn to a sharp increase 
at about the same temperature. This comparison implies 
that the relatively large values of protein reorganization 
energy obtained in our simulations at room temperature 
can be traced back to highly anharmonic motions of the 
protein matrix. 
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FIG. 8: Distribution of the population number of the ionized 
state of the special pair along the simulation trajectory at dif- 
ferent temperatures. The length of the simulation trajectory 
varies from 5 ns at 77 K to 15 ns at 300, 325, and 350 K. 



B. Polarizable special pair 

We need to emphasize here that our modeling of the 
polarizability of the special pair carries qualitative signif- 
icance only. In addition to the obvious limitations of the 
two-state model, the modeling of the temperature depen- 
dence of the special pair polarizability is not adequate. 
In our current simulations, the temperature dependence 
of the average population of the ionized, charge-transfer 
state of the special pair, ncT{T), originates from the tem- 
perature dependence of the diabatic diagonal energy gap 
fea l37p . This component of the two-state Hamiltonian in- 
creases with growing temperature, in general agreement 
with the idea that a polar environment should become 
effectively less polar with increasing temperature. There- 
fore, as is illustrated in FigureO the special pair becomes 
effectively more localized at higher temperatures because 
the average energy splitting between the two state grows 
with increasing temperature. The broad distribution of 
TT-CT is a signature of the strong vibronic coupling of the 
charge-transfer state»24 What effectively happens due to 
strong temperature dependence of the average popula- 
tion is that the polarizability of the special pair is about 
400 A3 at T = 300 K increasing up to 1800 at 77 
K. Given the experimental temperature variation of the 
absorption band of the special pair— and the results of 
Stark spectroscopy at 77 the former values appears 
to be more realistic than the latter. 

The increase of localization of the special pair in 
its neutral {Pm-Pl)* state results in a blue shift of 
the absorption spectrum in general agreement with 
experiment However, the slope of this temperature 
dependence derived from the data shown in Figure [5] ap- 
pears to be too large. In vibronic models of the tem- 
perature effect on the special pair absorption this ef- 
fect is modeled by temperature-dependent population 
of vibronic modes coupled to the dimer charge-transfer 
state.— This implies the temperature shift of the di- 
abatic diagonal energy gap. Since this property is de- 
termined by the protein/water electrostatic potential in 



our simulations, a possible way to off-set a too strong 
temperature dependence of absorption frequency is to in- 
troduce a temperature-dependent off-diagonal coupling J 
(Appendix [B|) Low frequency vibrations of the special 
dimer in the 90-160 cm~^ region^^ might contribute to 
that temperature dependence. It seems that the model 
needs to be modified to reproduce the temperature vari- 
ation of the absorption spectrum of the special pair. The 
current simulations in S2 protocol are therefore not ca- 
pable to properly address the issue of the temperature 
dependence of the rate. However, we still believe that 
our results provide valuable insights into how the param- 
eters of the model change once the polarizability is turned 
on. We therefore report the results of simulations here 
with the warning that the parameter magnitudes might 
be modified with the refinement of the model. We will 
also limit our analysis of the free energy surfaces of elec- 
tron transfer to 300 K at which the polarizability seems 
to be more realistic. What this value at room tempera- 
ture should be is not entirely clear since the Stark data 
were reported at 77 (see Appendix IB|) . 



C. Free energy surfaces 

A general solution for the non-ergodic free energy sur- 
face defined by eq[33]is still missing. The current calcu- 
lations and analysis of MD data are therefore limited to 
the phenomenological approach outlined in sec lll Cl where 
a step-wise frequency filter was introduced into the fre- 
quency linear response functions. Computer simulations 
and comparison to optical experiments in glass-forming 
liquids support this approach^ and one therefore can 
ask what would be the free energy surface G(A:et, X) on 
the time-scale of primary charge separation tet = ^et 
compared to the thermodynamic surface G{X). 

The thermodynamic free energy surface is of course 
not available to us since sampling is always an issue with 
simulations. However, leveling off of the reorganization 
energies on the 10-15 ns trajectory seen in Figure [5] al- 
lows us to hope that, except for the slowest modes re- 
sponsible for the conformational mobility of the protein, 
the phase space relevant to activating charge separation 
was adequately sampled. The free energy surfaces for 
non-polarizable (SI) simulations obtained from the 10 
ns trajectory for the initial (P-B^)* state and from the 
5 ns trajectory for the final (P"*"— B£)* state are shown 
in Figure m The results of polarizable (S2) simulations 
are collected in Figure [TUl Our simulations allow us to 
sample only the total interaction between the cofactors 
and the protein/water solvent (eq [5]) and therefore the 
gas-phase energy gap is missing from the overall energy 
gap X. This component of the energy gap was obtained 
from fitting the calculated rate constants at 300 K to 
the experimental data by Fleming et al^ and Wang et 
alM- (see below). The gas-phase gap obtained from the 
fit AES'^^ = 1.86 eV was used to horizontally shift G{X) 
obtained from simulations resulting in the average energy 
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FIG. 9: Free energy surfaces of primary charge separation 
obtained from MD simulations of reaction center with non- 
polarizable (constant charges, SI) special pair. The free en- 
ergy surfaces f3G{X) = — In P{X) have been obtained from 
the normahzed distributions of the total energy gap X from 
MD simulations (points) in the initial charge-transfer state, 
P-Bz, (10 ns trajectory), and the final state, P+-B£ (5 ns 
trajectory). The pair of curves marked with 0.21 are the non- 
ergodic free energy surfaces calculated by using f^^ = 0.21 
value of the non-ergodicity parameter following from the fit 
of the theoretical rate to experiment (T = 300 K) . The ver- 
tical separation of the initial and final free energy curves is 
-1-97 cm~^. The vertical separation of —450 cm~^ reported by 
Zinth and co-workers^^ is obtained when f^^ = 0.6 (marked 
in the plot) is used in the calculation of the final free en- 
ergy curve. The vertical separation between "equilibrium" 
curves is —1100 cm~^. The dash-dotted lines in the plot are 
fits to Marcus parabolas yielding equal reorganization ener- 
gies As ~ 1.6 eV consistent with direct calculations of second 
energy gap cumulants in Table [III The bold dashed line in- 
dicates the free energy obtained by solving the self-consistent 
non-ergodic equation for the rate (eg I28p by varying the aver- 
age energy gap (see the text) . The parameters are those used 
to calculate the charge-separation rate in Figure [TJ] 



gap of {AE) — 0.169 cV. This number, which is equal to 
the energetic separation of the free energy minimum from 
the point of activationless electron transfer, is consistent 
with the experimental value of 0.127-0.147 eV (from mu- 
tagenesis data) which separates the wild type reaction 
center from the top of the Marcus inverted parabolai^ 
Our result is also close to (AE) = 0.150 eV reported 
by Wang et alM: from fitting experimental data to the 
diffusion- kinetic model (see below). 

Long-trajectory simulations in SI protocol produce 
Marcus parabolas (dashed lines in Figure [5]) with the 
curvatures reproducing reorganization energies listed in 
Table [ITl Figure M also shows the non-ergodic parabo- 
las. Before explaining the calculation of those, we first 
need to comment on the experimental preparation of the 
initial state for charge separation. The initial state for 
primary charge separation is prepared by photocxcita- 
tion of the special pair which prior to that stays in the 
ground state for a time long compared to any relaxation 
time in the system. The ground state is thus character- 
ized by the equilibrium polarization P^q -I- Peq of which 
P' is the result of the inhomogeneous protein/ water en- 



FIG. 10: Free energy surfaces of primary charge separation 
obtained from MD simulations of the reaction center with 
polarizable (fiuctuating charges, S2) special pair. The upper 
curve is obtained from the simulation analysis with a 4 ps 
observation window, while the lower curve refers to the ob- 
servation window of 15 ns. The dashed line is the fit of the 15 
ns simulation data to the analytical Q-model with the fitting 
parameters: (AE) = 0.07 eV, Xs = 2.81 eV, and a = 0.45. 



vironment and Poq comes from the polarization of the 
environment by the special pair. When lifted to the 
excited state by the absorbed photon, the special pair 
changes its charge distribution and the polarization Poq 
can dynamically adjust to the new equilibrium polariza- 
tion P*q. We will assume that this change, P*q — Peq 
is insignificant compared to P^q -I- Poq on the reaction 
time-scale. We will therefore neglect the non-ergodicity 
correction in the Coulomb component of the shift assum- 
ing Ai?^(fcET) = AE'~" . This approximation results in 
the following non-ergodic free energy surface 



G(fcET,^) 



iX-{AE))' 

4A(fcET) 



In this equation. 



A(fcET) = A'"'^ + ./4(fcET)AC 



(45) 



(46) 



is the non-crgodic reorganization energy affected by the 
dynamical arrest of the Coulomb component of the sol- 
vent reorganization. The fit of the rate constant to ex- 
periment (see below) results in /,^o = 0.21, and the free 
energy surface obtained by using this non-crgodicity pa- 
rameter is shown by the solid line in Figure [5) 

There is a significant difference between the way the 
initial and final states for the first electron hop are cre- 
ated. The final state is characterized by an instanta- 
neously created dipole moment of the charge-separated 
state and so the corresponding Stokes shift requires non- 
ergodic correction with the following result for the final 
free energy surface 



G'{X) = 



(X - (AE) + fLAX^t) 
4A' 



(47) 



In this equation, AX^t is the total Stokes shift between 
the minima of two parabolas achieved on long simulation 
trajectories. The non-ergodicity parameter f^^ and the 
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reorganization energy A' depend on both the life-time of 
the charge-separated state and the corresponding Stokes 
shift dynamics. FinaUy, the vertical energetic separation 
between the parabolas' minima AGne (not the reaction 
free energy) follows from the condition G'(fcET, 0) = G'(0) 
once all other parameters are known. 

We currently do not have sufficient data to calcu- 
late the non-ergodic parameters in eg 1471 (which require, 
among other things, the free energy surface correspond- 
ing to the electron located at H^) and so will limit our 
arguments to qualitative considerations only. The equi- 
librium free energy surfaces obtained from long simula- 
tion trajectories are vertically shifted by AG = —1100 
cm~^. This number is consistent with experimental data 
from recombination ratesi^J^ which have put the lowest 
limit for AG at ~ —2000 cm~^. This later value might be 
overestimated since it was measured on the 100-/is life- 
time of the triplet state of the special pair. Delayed flu- 
orescent measurements^^ sampling the system on the 20 
ns time-scale and photovoltage measurements at the 15 
ns time-scale^SS both comparable to the length of sim- 
ulations, show somewhat smaller gaps, AGno — —1370 
cm~^ and —1180 cm~^, respectively. The latter data 
refer, however, to the Rhodo spirillum rubrum reaction 
center. 

We need an assignment of /^o in eq[47]to produce the 
non-ergodic free energy surface of the charge-separated 
state. If we use /^^ = f^^ = 0.21 from the analysis of 
the primary charge separation rate, we get essentially no 
vertical shift of the two parabolas, AGne = 97 cm~^. 
The vertical shift of —450 cm~^ reported by Zinth and 
co-workers2^ is obtained when /^^ = 0.6 is used in ealTTl 
This latter value of the vertical displacement of parabo- 
las minima, measured on the picosecond time-scale, com- 
pares well with the estimate by Holzwarth and Miillerj^S 
—331 cm~^, also done on the picosecond scale. We will 
postpone a more detailed analysis of the energetics of 
subsequent electron hops to a future publication, while 
the current analysis is aimed to show that overall our 
results do not contradict the key experimental observa- 
tions reported in the literature. We only note here in 
passing that, similar to our previous simulations of hy- 
drated plastocyaniufiSi our present simulations show a 
clear separation between the Stokes shift AXgt and twice 
the solvent reorganization energy, 2Xs (also see ref |65|). 
We will address this problem in more detail elsewherc^^^ 

We need to caution here against a too literal under- 
standing of the non-ergodic free energy surfaces of elec- 
tron transfer. Under ergodic conditions, the free energy 
surface can be sampled by changing the average energy 
gap by, for instance, optical spectroscopy. The result can 
then be directly applied to the Frank-Condon factor of 
the reaction yielding the reaction energy gap law. In the 
case of non-ergodic reactions, this direct application of 
the free energy surface obtained at a given observation 
window is prohibited since the spectrum of fluctuations 
changes at each rate constant achieved by horizontally 
sliding the free energy surface and thus sampling the av- 



erage gap. In order to illustrate that, we have plotted 
in Figure [9] the free energy surface obtained by changing 
the average energy gap in the self-consistent non-ergodic 
equation for the rate constant (eg 1^5)). The result is a 
funnel-like surface, which we also previously obtained in 
a study of ergodicity breaking in liquid crystalsi^S^ The 
such obtained curve transforms from the narrow free en- 
ergy surface at a high reaction rate to the thermodynamic 
surface when the barrier for the reaction increases and 
the rate slows down. 

As is clear from the broader distribution of energy gaps 
for the polarizablc special pair (Figure and from Table 
|TT] where specific values of the reorganization energies are 
listed, the free energy surfaces G{X) are quite different 
for a polarizable and non-polarizable special pair. As a 
matter of fact, not only curvatures (reorganization en- 
ergies) are different in two cases, but also the shape of 
the free energy surface changes from a Marcus parabola 
in the former case to a significantly asymmetric shape in 
the latter (Figure [TU)) . This result is consistent with the 
predictions of the Q-model of electron transfer in polariz- 
able donor-acceptor complexes and, in fact, the simulated 
curve is well fitted to eqdl] (dashed line in Figure [TU)) . 
Note that the reorganization energy obtained from the fit 
is close to the result of direct calculation from the second 
cumulant (eq[T5l Table [Til . 

The increase in the reorganization energy in the case of 
polarizablc P comes from fluctuations of the amount of 
charge transfer between the covalent and ionized states 
of P (Figure [TT|) . It is clearly seen from Figure [11] that 
energy gap fluctuations in excess to those existing for 
non-polarizable P trace the fluctuations of ncT- Most 
of the excess reorganization energy comes from the pro- 
tein. The reorganization energy from water actually gets 
smaller when polarizability is introduced, but the protein 
reorganization energy is increased by a factor of four. 

We do not currently have an established theoretical al- 
gorithm of how to calculate the non-parabolic free energy 
surfaces of electron transfer involving highly polarizable 
donor-acceptor states when ergodicity is broken. In the 
absence of a theoretical formalism, we have turned to 
simulations. Figure [TUl shows the free energy surface pro- 
duced from simulations by sliding the observation win- 
dow of the length 4 ps along the trajectory and then 
averaging all the histograms produced from each window 
after shifting them to a common probability maximum. 
The normalized distribution produced in this way is then 
used to plot the non-ergodic free energy curve shown in 
Figure [TUl In contrast to distributions obtained with the 
non-polarizable simulation protocol, the non-ergodic sur- 
face turns out to be non-parabolic. We do not presently 
have a good explanation of this observation. 



D. Charge-transfer rates 

The decay of the population P{t) of the photoexcited 
special pair is known to be non-exponential4^ii2iii£^ Re- 



15 




3 -1 



Polarizable 

(b) 

Non-polarizable 




2000 



4000 6000 
t/ps 



8000 



10000 



FIG. 11: The trajectory of the population of the charge trans- 
fer state of P (a) and the trajectory of energy gap fluctuations 
(b) for non-polarizable (black line) and polarizable (gray line) 
special pair. 



cent explanations of this observation^ii2£ have cast the 
problem in terms of the Fokker-Planck kinetics with a 
Golden Rule reaction sink, similar to formalisms devel- 
oped in the past by Agmon and HopfieldiSi and Sumi 
and Marcus4^ This theoretical algorithm offers the fol- 
lowing physical picture. At the initial time t = 0, a 
laser flash lifts the equilibrium population Pcq{X) of the 
ground P to the excited state P* (dashed line in the left 
panel in Figure [T2)) . At this moment, the state P* is 
fully occupied, P(0) — 1. The system can decay to the 
charge-separated state with the frequency oje (eqH)) at the 
activated state X = thus depleting P{t) and changing 
Pcq{X) to P{X,t) (dash-dotted line in the left panel in 
Figure [T^. At the initial time, P{X,t) ~ Pcq{X), and 
the decay is determined by the equilibrium rate /cet given 
by eqISl However, as the population of the activated state 
X = depletes from that given by the Boltzmann dis- 
tribution, the continuation of the reaction requires a dif- 
fusional supply of the population to the activated state. 
The result is a slower population decay and effectively 
multiexponential kinetics. Given that the activation bar- 
rier is small for primary charge separation (Figure [H]), 
the diffusional regime kicks in at the early stage of the 
reaction leading to observable deviations from monoex- 
ponentiality. 

Two complications need to be recognized in applying 
this type of diffusion-reaction kinetics to the problem of 
primary charge separation. The first complication, well- 
recognized in studies of the dynamic solvent effect on 
electron transfer in small molecules j^iiS^ is related to the 
fact that the Stokes shift dynamics are non-exponential, 
in particular in its initial Gaussian stage. The common 
approach to the problem, going back to the Sumi-Marcus 
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FIG. 12: Left panel: Photoexcitation of the special pair lift- 
ing the equilibrium distribution (dashed line) from the ground 
state to the electronically excited state. This excitation starts 
the decay of the population through the Gaussian sink k{X), 
along with the one-dimensional diffusion given by the Fokker- 
Planck operator LP{x,t). The dash-dotted line indicates de- 
pletion of the population at the side of the sink resulting in a 
slowing down of the population relaxation and in overall non- 
exponential kinetics. Right panel: Population decays of mu- 
tants of the reaction center of Rhodobacter sphaeroides taken 
from refill (points) and fits of P{t) to the diffusion-reaction 
model (solid lines). The legend in the right panel specifies mu- 
tants altering the local electrostatic potential at the special 
pair.— The dashed line marked /cet shows the initial popula- 
tion decay with the electron transfer rate constant according 
to eq[S21 



formalism)^ is to split the overall energy gap X into a 
fast, Xf, and slow, a;, components, X = Xf + x. The 
evolution operator along the reaction coordinate X is 
then averaged over the equilibrium distribution of the 
fast component, resulting in a diffusion-reaction equation 
for the population along the slow reaction coordinate x 
only: 

dP{x,t)/dt= [L{kET,x) - k{x)]P{x,t) (48) 
In this equation, L{kET,x) is a diffusional operator 



L{kET,x) = D{kET) 



d_ 

dx 



d_ 

dx 



dG{kET, x) 
dx 



(49) 



describing the Fokker-Planck dynamics in the poten- 
tial given by the electron-transfer free energy surface. 
For multi-exponential decay, a time-dependent diffu- 
sion constant can be used for the harmonic potential 
G(fcET, a;),"i2£ while an effective relaxation time Tes needs 
to be defined for a general potential. Following HyneSfiiS 
this relaxation time is defined here in terms of a weighted 
sum of the corresponding rates of exponential relaxation. 
For a bi-exponential long-time tail in eg 1211 one gets 



'off 



(Afrr^+A^r,-i)/(Af + A^) 



(50) 



The diffusion constant in eg then becomes D{kET) = 
2k^TX'~'{kET)/TeS- Finally the rate constant k{x) in eq 
25] is the Golden Rule rate averaged over the equilibrium 
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distribution of the fast relaxation component 



k{x) = ujeJXs/{>^G + ^'"'') cxp 



4(Ag + Ai°d) 



(51) 

where Xq is the fast Gaussian component of decay in eq 
[21] A'"'^ is the induction reorganization energy, and (Ai?) 
is given by eq\5\ 

Most studies applying this formalism in the past have 
assumed that the overall rate of diffusional reaction, i.e. 
the rate of arriving at the transition state X = from 
the bottom of the potential well, is much smaller than 
the relaxation rate of any nuclear mode coupled to elec- 
tron transfer. This is obviously not true in our case, 
and non-ergodicity corrections, already introduced into 
eqs l48l - [5n are required. These corrections come in the 
form of the free energy surface G{k-ET,x) depending on 
the rate /cet (eq[34]), as well as the diffusion coefficient 
D{kET) depending on the non-ergodic reorganization en- 
ergy A*-^(fcET)- Therefore, any solution of the dynamic 
diffusion-reaction equation should produce a closure for 
^ET and then solved by repeated iterations.^^^ Equation 
25] was solved employing the generalized moment expan- 
sion of Nadler and Marcus^ to produce /cet as the initial 
population decay (right panel in Figure [9]) 



dlnP{t) 



dt 



(52) 



where P{t) = J P{x,t)dx. This condition establishes the 
closure for the self-consistent calculation of ^et by re- 
peated solutions of eq |48l The free energy surface is 
then obtained by a horizontal shift of eq|45l G'(fcET, x) = 

xV[4A(fcET)]. 

The approach outlined here results in a good agree- 
ment with experimental population decays for a number 
of mutants reported by Wang et al^ (Figure [T^ with 
the input parameters produced by SI simulation proto- 
col. Also notice that the rate constant in the sink term in 
easl48land[511is purely classical and does not incorporate 
quantum vibrations. For reactions in the inverted region, 
quantum Franck-Condon vibrational overlaps provide ad- 
ditional vibronic channels for electronic transitions^^ 
Primary charge separation appears to operate in the nor- 
mal regionS^ (Figure [5]) and quantum vibrations can be 
omitted. Notice, however, that classical phonon modes 
have been included into the fast Gaussian and induction 
components of the reorganization energy. 

The Stokes shift correlation function necessary to cal- 
culate the non-ergodic reorganization energy from eg 1271 
at each iteration step in eq 05] was taken from MD sim- 
ulations of the reaction complex (Figure [T3|) . Several 
important observations follow from examining Figure [T51 
The ballistic component of the decay, arising from bal- 
listic motions of water and quasi-lattice vibrations of the 
protein matrix, is significantly diminishedi^ compared 
to the Stokes shift dynamics of small chromophores in 
waterji^ Indeed, the sum of the Gaussian component of 
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FIG. 13: Stokes shift correlation function of primary charge 
separation obtained from 5-15 ns MD trajectories in SI and 
S2 simulation protocols. S{t) at different temperatures with 
non-polarizable special pair are shown in (a). In (b), the over- 
all Stokes shift correlation function at 300 K (marked as X) is 
compared to its components from Coulomb and induction in- 
teractions, along with the Coulomb/induction cross term (eq 
I17[l . In (c), the same separation into components of S{t) is 
shown for simulations with the polarizable special pair at 300 
K. 



the Coulomb reorganization energy and the induction re- 
organization energy, both responsible for the fast decay, 
is below 20% of the overall solvent reorganization energy 
As. This fact is critical for the analysis of non-ergodic 
free energy surfaces of electron transfer as the fast relax- 
ation component is essentially the only part of nuclear 
reorganization which is not dynamically arrested on the 
short time-scale of charge separation (see below). 

The exponential tail of the Stokes shift decay becomes 
slower with cooling, as expected. The effective relaxation 
time roff(T) can be calculated from the fitted exponential 
relaxation times according to eqs [22l and [50l and fits well 
by an Arrhenius function (200 < T < 400 K) with the 
activation energy of Er = 1212 K. This activation barrier 
is close to the value ~ 1060 K reported for the long-time 
tail of the fluorescence decay of an optical probe bound 
to a proteiufii^ This activation barrier was assigned to 
local segmental motions of the protein coupled to the 
hydration layer. The long tail in the Stokes shift corre- 
lation function observed here is, however, shorter than 
that reported in ref Ill3l and is in fact close to the slow 
protein-water dynamics with the characteristic time of 



^90 ps recently reported from Stokes shift data in ref 

The combination of the Arrhenius temperature depen- 
dence with the low activation energy points to the link 
between exponential Stokes shift relaxation and /3 relax- 
ation of the protein/water system. Previous measure- 
ments of a relaxation in hydrated proteins have consis- 
tently shown much larger effective activation barriers of 
the order 6000-9000 K^^^ in addition to the breaking 
of the Arrhenius law in a broad temperature range. We 
therefore conclude that primary charge separation is cou- 
pled to two nuclear modes: Gaussian ballistic/phonon 
motions and exponential /3 relaxation. The relaxation 
time of the former turns out be close to 0.1 ps^ and is 
essentially independent of temperature. We note in this 
regard that anharmonic protein displacements shown in 
Figure[7]are also linked to /3 relaxationi^ The decoupling 
of the Stokes shift dynamics of the primary charge separa- 
tion from a relaxation is distinct from the situation com- 
monly seen for solvation dynamics of small solutesii and, 
among other things, implies that dielectric a- relaxation 
data, routinely used to calculate solvation dynamics of 
small chromophoreS)^ have little to do with the dynam- 
ics of primary charge separation. 

Stokes shift dynamics allow us to calculate the non- 
ergodic reorganization energies entering the reaction 
rates and population decays. Table [TTIl reveals yet an- 
other important mechanistic aspect. It shows that the 
reorganization energy of water is significantly cut off by 
the dynamical arrest. Reorganization of fast, anharmonic 
quasi-lattice vibrations of the protein emerges from the 
water dominance in the thermodynamic limit, acting as 
the leading mode driving electronic transitions on the pi- 
coseconds scale. 

We now turn back to the calculation of the rates of pri- 
mary charge separation. Two types of laboratory exper- 
iments are most relevant to our discussion. The first are 
the measurements by Fleming and co-workers^ of charge 
separation rates in a broad range of temperatures be- 
tween helium 5 K and room temperature, 300 K. The ex- 
perimental observation, which has puzzled theorists ever 
since, is a very gentle decay of the electron transfer rate 
over the whole temperature range (open points in Fig- 
ure I14p . This result is apparently inconsistent with any 
conceivable temperature dependence of equilibrium nu- 
clear solvation energies, even if activationless electronic 
transition is realized at some intermediate temperature. 
The second set of experimental results, reported by Allen 
and Woodbury and co- workers, ^^^iSS provides population 
decays of P* in a carefully engineered set of mutants al- 
tering the electrostatic potential at the location of the 
special pair. A surprising result of these experiments was 
the realization that the wild-type reaction center falls off 
the top of the energy gap law into the normal region of 
electron transfer i2S 

Our current calculations, based on the input from MD 
simulations and the concept of solvation non-ergodicity, 
are capable of reproducing the experimental decay curves 
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FIG. 14: Temperature dependence of the rate of primary 
charge separation from experiments by Fleming et aJi (points) 
and from calculations of the rate using the gas-phase energy 
^£;gas ^ -j^ gg gY y ^ ^-j^ g ^j^-i adjusted to reproduce 
the rate at T = 300 K (solid line). The dotted line repre- 
sents the fit of experimental data to an empirical equation 
suggested in ref [J. 



P{t) for the whole set of mutants studied by Allen and 
Woodbury (Figure fT2|) . In the fit, the electron transfer 
matrix element was obtained from the 300 K rate of the 
wild-type reaction center and the inhomogeneous compo- 
nent of the average energy gap AE^^^ was varied among 
the mutants (electrostatic mutationals). The fitted vari- 
ation of AE^^^^ is consistent with changes in the midpoint 
electrochemical potential upon the mutation (Appendix 
\G\\ . In order to further test the consistency of these re- 
sults with the experimental database, one needs to prove 
that the experimental rates at different temperatures^ 
can be obtained with the set of parameters used to fit 
the mutagenesis data. These results are shown in Figure 
[14] with the details of calculations given in Appendix ICl 

Proper account of the temperature variation of the pa- 
rameters entering the activation barrier is important in 
reproducing the observed rates. The main factor here 
is the temperature dependence of the induction shift, 
which is well converged in our simulations and slopes 
positively with increasing temperature (Table |T| . Unfor- 
tunately, the accuracy of the current simulations does 
not allow us to address the temperature dependence of 
the Coulomb components of the energy shift and reor- 
ganization energy since their changes in the interval of 
temperatures studied are within the uncertainties of nu- 
merical simulations. Our previous experience with an- 
other photosynthetic protein, plastocyanin, suggests that 
the length of the simulated trajectories needs to be ex- 
tended up to at least 20 ns for a reliable estimate of 
the temperature slope,-^^ which is beyond our current 
computational capabilities. Therefore, in order to as- 
sign realistic slopes to the Coulomb components of the 
free energy barrier, we used our previous observatio n^^i^^ 
that the results of many calculations and experiments 
on small solutes in polar solvents give the logarithmic 
slope of the Coulomb reorganization energy in the range 
AA'^'/A'^ ~ -(2-3)xlO-3AT. The Coulomb reorganiza- 
tion energy was then assigned the slope of A In /AT = 
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TABLE III: Solvent reorganization energies (eV) and their water and protein components from MD simulations. The reorga- 
nization energies depending on the reaction rate were obtained from Stokes shift dynamics according to eg 1271 /cet refers to 
the charge separation rate at the corresponding temperature. 



Protocol T/K 


A. 


Aa(fcET) 


Aw 


Aw(fcET) 


Aprot 


Aprot(fcET) 


SI 


77 


0.421 


0.266 


0.191 


0.072 


0.245 


0.210 




200 


0.903 


0.257 


0.756 


0.106 


0.251 


0.201 




250 


1.955 


0.261 


1.639 


0.126 


0.419 


0.235 




300 


1.598 


0.454" 


1.124 


0.124 


0.466 


0.260 




350 


2.239 


0.657 


1.246 


0.202 


1.407 


0.611 




400 


1.335 


0.639 


0.866 


0.259 


0.797 


0.451 


S2 


300 


2.513 


1.276 


1.454 


0.188 


1.385 


0.970 



"Direct fits of the experimental population decays to the Sumi- 
Marcus model considering the average energy gap and the reor- 
ganization energy as fitting parameters gave As = 0.350 eV and 
(AE) = 0.150 eVi^ 
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FIG. 15: Charge separation rate vs the variation of the av- 
erage donor-acceptor energy gap produced by mutagenetic 
substitution (points^). The lines are obtained by horizontal 
shifts of the non-ergodic parabolas of the initial charge sepa- 
ration state obtained in SI protocol (solid line, Figure[9ll and 
in S2 protocol (dashed line, Figure [TU]). In experiment, mu- 
tagenetic substitution varies the inhomogeneous part of the 
Coulomb component of the average vertical gap and there- 
fore that parameter marks the horizontal axis, AAE^-^^^ = 
corresponds to the wild-type reaction center. 

— 1.3x 10^"^ K^^ and, based on its relative magnitude, the 
Coulomb component of the average energy gap was given 
the slope of A(ln A£;'^)/Ar = 5.2 x lO"" K'^ (see Ta- 
ble IIVI in Appendix [C]) . These assignments do not affect 
our results much since a close fit can also be obtained by 
assuming these these two parameters are temperature- 
independent. We finally note that the problem of the 
temperature dependence of the reaction parameters, in 
particular the driving force, is not free of controversy. 
Opposite signs of reaction entropy have been obtained 
for different charge-transfer reaction a^^i^^^ in the reac- 
tion center and temperature-independent parameters are 
routinely used in the analysis F^i^iii^ 

The non-ergodic free energy surface of electron trans- 
fer (narrower surface in Figure [TU|) obtained from simula- 
tions with a polarizable special pair can also be used to fit 
the experimental reaction rate at 300 K. This fit results 
in the gas-phase gap of AE^^^ = 1.57 eV used in Figure 
[TUlto plot the free-energy surfaces. This value of the gas- 



phase gap yields the average energy gap of (AE') = 0.150 
eV, consistent with the experimental evidence^^ and pre- 
vious fits of the rates by Wang et alM The non-ergodic 
reorganization energy obtained from fitting the curvature 
at the minimum of the G(fcET,^) curve turns out to be 
0.39 eV, close to the value of 0.46 eV reported for SI sim- 
ulations in Table IIIII and the value of 0.35 eV reported 
by Wang et al^ The non-ergodic free energy curves from 
Figures [9land[T0l arc used to construct the energy gap law 
of electron transfer plotted against the variation of the 
inhomogeneous electrostatic potential of the protein, as 
was done in mutagenesis experimentsi^ We find that po- 
larizable and non-polarizable simulations result in close 
shapes of the energy-gap law. 



IV. DISCUSSION 

A. Mechanism of electron transfer activation 

The extensive MD simulations of the bacterial reac- 
tion center combined with formal modeling have allowed 
us to look closely at the nuclear modes driving electronic 
transitions and their energetic balance in the reaction ac- 
tivation barrier. Several qualitative results have emerged 
from our analysis. From the viewpoint of the relative par- 
ticipation of different types of interaction potentials, we 
have shown that induction and Coulomb forces give com- 
parable contributions to the average energy gap, while 
Coulomb interactions tend to dominate the reorganiza- 
tion energy of electron transfer. A significant finding 
of this study is the realization that, on the nanosecond 
time-scale achievable by numerical simulations, the reor- 
ganization energies and shifts are quite significant, much 
larger than had been anticipated so far. The under- 
standing that most of this nuclear solvation is dynam- 
ically frozen on the time-scale of the reaction then be- 
came critical for the quantitative description of the ob- 
servable rates. While water dominates the reorganiza- 
tion energy on the nanosecond time-scale, most of this 
solvation freezes on the picosecond reaction time-scale. 



19 




20 40 60 80 100 

t /ps 



FIG. 16: Total Stokes shift correlation function of primary 
charge separation compared to the Stokes shift correlation of 
reaction center tryptophans (dashed line) and to the normal- 
ized self-correlation function of the fluctuations in the donor- 
acceptor distance RuA{t) between the special pair and bacte- 
riochlorophyll (dash-dotted line). The lower solid line shows 
the correlation function of tryptophan absorbance band taken 
from ref |4^. 

and protein vibrations emerge as the main nuclear mode 
driving electronic transition. Still, there is a noticeable 
component of water reorganization, originating from the 
ballistic Gaussian decay of the Stokes shift correlation 
function, left even on the picosecond time-scale (Table 

imi. 

Once the protein is identified as the major heat reser- 
voir operating on the picosecond time-scale of the re- 
action, one can try to identify a particular mode most 
strongly coupled to the transferred electron and driving 
the electronic transition. Several answers to this ques- 
tion have been proposed in the past. Wang et al^ sug- 
gested to use transient changes in tryptophan absorbance 
at 280 nm to monitor the electron transfer kinetics. In 
this approach, photoexcited tryptophan serves as a time- 
resolved probe of the ultrafast nuclear rearrangement 
of the protein matrix with the hope that the dynamics 
recordcrcd by spectroscopy will match the Stokes shift 
dynamics unreachable by spectroscopic techniques. Since 
both types of information arc available from our simula- 
tions, we have tested this hypothesis here. 

Figure [12] compares the Stokes shift dynamics of pri- 
mary charge separation with the Stokes shift dynamics of 
tryptophan averaged over all tryptophans in the reaction 
center protein. Although these two match each other 
reasonably well, the experimental traced shown in the 
same plot is quite different having, in particular, a much 
faster initial decayJ^ It turns out that this experimental 
trace matches quite well the autocorrelation function of 
the donor-acceptor distance between the special pair and 
the accessory bactcriochlorophyll cofactor B l (Figure [51 
also see Figures SI and S2). The decay of this function is 
also caused by protein vibrations suggesting that the ex- 
perimental observation traces one of the long-wavelength 
vibrational modes responsible for large-scale protein mo- 
tions, but not necessarily the modes contributing pri- 
marily to the Stokes shift dynamics of primary charge 
separation. 

In fact, following an early proposal by Gehlen et 



al^ Chaudhury et alJ^ recently suggested that donor- 
acceptor vibrations represent the mode activating elec- 
tronic transitions. Our current results do not support 
this hypothesis. The dynamics of the donor-acceptor 
vibrations are different from the Stokes shift dynamics. 
In addition, both the self-correlation function of donor- 
acceptor distances and the experimental trace of Wang et 
al^ produce too large an amplitude of the initial decay 
which would make a larger portion of nuclear solvation 
unfrozen on the time-scale of the reaction, thus invalidat- 
ing the analysis of the reaction rates (see below). On the 
experimental side, related evidence shows the low sensi- 
tivity of charge-recombination rates to high pressure (up 
to 345 MPa) which caused about 16% of volume change 
of the samplcii^i 

What has not been considered so far in the long history 
of modeling primary charge separation is the possibility 
that a high polarizability of the photoexcited special pair 
can significantly modify the energetics of the reaction. 
Our simulations here are the first attempt to understand 
the possible consequences of the gigantic polarizability of 
the special pair for the charge-transfer energetics and ki- 
netics. The polarizability of the special pair was modeled 
here by the two-state mode l^^i^^^ with the dynamic ad- 
justment of the population between charge-transfer and 
neutral states of the primary pair along the simulation 
trajectory. The two-state model has its obvious limita- 
tions and a possibility of a broader spectrum of electronic 
states^^*^ can be considered in the future, along with 
improved modeling of the temperature variation of the 
absorption spectrum of the special pair. Nevertheless, 
the present simulations give the first insights into what 
sort of changes to the energetics should be anticipated 
when the polarizability has been taken into account. 

What we have found here is consistent with previous 
studies of the role of polarizability in the energetics of 
electron transfer The free energy surface of the ini- 
tial electron-transfer state involving a polarizable spe- 
cial pair is significantly distorted compared to Marcus 
parabolas which we obtained in the simulation protocol 
with the non-polarizablc primary pair (c/. Figures [5| and 
nop . The reorganization energy, obtained as the variance 
of the energy-gap fluctuations, is significantly enhanced 
compared to the case of non-polarizable simulations, also 
in agreement with the previous studies The free en- 

ergy curve could be fitted with the analytical equations 
of the Q-model, which introduces a non-parabolicity pa- 
rameter in addition to the two-parameter description of 
the Marcus model. 

Although the free energy curve from the 15 ns tra- 
jectory shown in Figure [TUj is perhaps the most asym- 
metric electron transfer free energy surface ever reported 
from numerical simulations^^ most of this asymmetry is 
washed out by the dynamical arrest of nuclear solvation 
on the reaction time-scale. The free energy surface nar- 
rows down and approaches the Marcus parabola on the 
4 ps observation window (Figure fTU]) . In fact, the rate 
of charge separation can be equally well described by ei- 
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thcr polarizable or non-polarizable simulation data with 
a close range of parameters and a close match between 
the resulting energy gap laws (Figurc fTS)) . It appears that 
what charge separation probes on the picosecond obser- 
vation window is a stripped surrogate of the rich dynam- 
ics and thermodynamics of the protein/water electrostat- 
ics on the time-scale of thermodynamic observables. Na- 
ture has therefore played with dynamical time-scales to 
reduce these complexities to a near-resonance electron 
transfer driven by ballistic phonon motions. 



B. Rates of primary charge separation 

The ideas of non-ergodic nuclear solvation advocated 
here were tested for consistency with experimental ob- 
servations by calculating the rates of charge separa- 
tion as a function of tcmperaturCfi and the kinet- 
ics of the population decay depending on mutagenetic 
substitution!^ Since the simulation protocol involving 
the non-polarizable special pair had produced Marcus 
parabolas for the free energy surfaces, we were able to 
use the diffusion-kinetic model advanced by Sumi and 
Marcus^ in order to calculate the population decay. The 
main question we were asking here is whether the use 
of the same set of fitting parameters (electron-transfer 
matrix element V and the gas-phase energy gap AE^'^^) 
would provide us with a consistent description of both 
sets of experiments. We obtained a positive result here 
(Figures [n] and [HI). 

A fit of the population decay to the Sumi-Marcus 
diffusion-kinetic model was also presented by Wang 
et al^ In their analysis, the diffusion coefficient was 
taken to be time-dependent in order to reflect the non- 
Markovian character of the relaxations^ Both the aver- 
age energy gap and the reorganization energy were con- 
sidered as fitting parameters. We found that the use of 
the time-dependent diffusion coefficient is not a necessity 
since the same data can be reproduced with an effec- 
tive diffusion coefficient extracted from Stokes shift cor- 
relation functions. What distinguishes our analysis from 
their's is that the solvent-induced shift of the average en- 
ergy gap and the reorganization energy are fixed by MD 
simulations, instead of used as fitting parameters. The 
gas-phase gap and the electron transfer matrix element 
were fitted to the rate at 300 K, but then, these parame- 
ters allowed us to reproduce Fleming's data. With these 
restraints on the parameters' magnitudes, there is very 
little room for adjusting the two parameters. We also 
note that the results of the calculations, and of our sim- 
ulations of the charge-separated state P+-B^, are con- 
sistent with the current state of experimental evidence 
regarding the energetics of primary charge separation, as 
we have discussed in sec IIIICI 

The negative slope of the charge separation rate with 
increasing temperature has puzzled theorists for two 
decades, and has mostly been approached by considering 
a temperature-dependent population of phonon modes 



coupled to electron transfer Although our simulations 
and conclusions are for the most part limited to high 
temperatures greater than 200 K, explaining the nega- 
tive temperature slope of the rate in this region does not 
require vibronic coupling models. We found the reaction 
rate to follow the temperature variation of the induction 
component of the average energy gap which itself be- 
comes less negative with increasing temperature because 
of the protein expansion. 

We have confirmed the observation made by Haffa et 
al^ that primary charge separation falls into the nor- 
mal region of electron transfer (Figures [9] and fTS]) . This 
result was considered incompatible with the weak tem- 
perature dependence of the rate, and a vibrational heat- 
ing mechanism^^ii^ii^ was suggested in order to explain 
the positive (Ai?). Our current calculations suggest that 
Fleming's data can be reconciled with (AS) ~ 0.15 eV 
for the wild-type reaction center without assuming vibra- 
tional heating once the temperature dependence of (Ai?) 
is taken into account. The main component of (Ai?) 
responsible for its temperature dependence is the shift 
by electronically instantaneous induction forces which 
do not get dynamically frozen, but rather change due 
to a temperature-affected alteration of protein's refrac- 
tive index (eq [33]) . Our simulations also suggest that 
the wild-type reaction center is driven even further from 
the optimum activationless energetics when temperature 
increases above the room temperature and that the op- 
timum activationless configuration is reached at around 
200 K (Figure S3). We refrain from speculations on evo- 
lutionary implementations of this result. 



C. Broader insights 

Electron transfer connects cofactors in energetic redox 
chains in biology. Three parameters are generally be- 
lieved to have the main impact on the kinetics of elec- 
tron hops: the redox potential, the probability of tun- 
neling, and the reorganization energy. The first one 
is relatively well understood, and in many cases, ac- 
cessible to measurements. The distance decay of elec- 
tron tunneling has attracted significant attention of the 
thcoreticali22ii^ and experimental^ communities in re- 
cent decades. Although the importance of specific path- 
ways in the polypeptide structure vs the generic tun- 
neling decay specified by the height of the potential 
barrier is still actively discussed )^^^'^'^°'^'^^i^'^^ there is a 
general consensus about the magnitudes of matrix ele- 
ments involved and the distance decay of the tunneling 
probabilityJ^Sii^ 

The last component of the biological electron transfer 
picture, the reorganization energy, is probably least un- 
derstood. Although the reorganization energy is the hall- 
mark of the classical Marcus theory of electron transfer^ 
not much is known about both its value and the micro- 
scopic modes responsible for reorganization in proteinic 
and DNAi^i electron transfer. For proteins, the exper- 
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imental evidence mostly comes from kinetic measure- 
ments of ruthenium-modified proteins introduced into 
the field by Gray and co-workers and some recent 
reports from computer simulations Notice that 
computer simulations reported in the past were mostly 
limited to either very short trajcctoric o^^i^°i^^^ or es- 
timates of the reorganization energy from the Stokes 
shift ^li^ which does not necessarily provide the correct 
value of the reorganization energy defined through the 
variance of the energy gapf^^i The uncertainties of reor- 
ganization energy values have led Button and co-workers 
to suggesti^S a generic value of 0.7 eV for electronic tran- 
sitions between cofactors not exposed to water with the 
provision that smaller values might be required for photo- 
synthetic electron transfer. The fits of the photosynthetic 
rates have been attempted many times and extremely 
low values of the reorganization energies (as low as 0.1 
eV— lii^), completely unthinkable in light of our present 
simulations, have been reported in the literature. 

Our present work gives a different perspective to the 
problem of the activation barrier of electron transfer in 
proteins. We claim that the range of reorganization ener- 
gies fundamentally attainable in protein electron trans- 
fer is very broad given that the overall reorganization 
energy attained in our present and previousi^i simula- 
tions is much higher than it was previously anticipated 
(~ 1.6 eV for SI protocol and ~ 2.5 eV for S2 proto- 
col). The question of assigning the reorganization energy 
thus turns not into its "generic" value, but into the ques- 
tion of finding the protein/solvent reorganization energy 
reachable on a given time-scale of the reaction, when a 
certain portion of nuclear degrees of freedom is dynami- 
cally frozen. 

We could not identify any specific solvent and/or pro- 
tein modes that drive electron transfer. Instead, the en- 
ergetics of electronic transitions appear to be driven by 
some generic set of ballistic modes which would prob- 
ably characterize any heterogeneous solvent made of a 
rigid core (protein) surrounded by a molecular polar sol- 
vent (water). It also seems true that achieving both 
the reaction rate of primary charge separation and its 
low temperature dependence allows some, although not 
large, fiexibility in the driving force (~ 0.3 eV between 
photosynthetic bacteriaF^). Where the specific design of 
the reaction center appears to matter is in providing a 
sufficient tunneling rate between closely separated cofac- 
tors. This part of the design turns out to be very essen- 
tial since the fast rate allows the natural photosynthesis 
to dynamically freeze nuclear solvation, and to achieve 
low values of the reorganization parameters character- 
ized by weak temperature dependence (ballistic motions 
and local density fluctuations). It might therefore turn 
out that "Darwin at the molecular scale"^^ operates not 
that much with redox potentials but, to a greater extent, 
with relaxation time-scales. 

Supporting Information Available: Atomic 
charges of the bacteriochlorophyll cofactors from DFT 
calculations and Figures S1-S3. This material is available 



free of charge via the Internet at http : //pubs . acs . org. 
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APPENDIX A: SIMULATION PROTOCOL 

Amber 8.0^ was used for all MD simulations and mini- 
mizations. The initial configuration of the reaction center 
complex was taken from a crystal structure of the purple 
bacterium Rhodobacter sphaeroides^ The force fields of 
bacteriochlorophylls, pheophytins, ubiquinones, and iron 
center were taken from Marchi and coworkersi^ The pro- 
tocol for the creation of solvated micelle was taken from 
ref m with the slight variations described below. 

First, it should be mentioned that the reaction center 
was built without the carotenoid cofactor, since it was 
deemed unnecessary for the photosynthesis functionJ^ 
The system was initially setup by protonating all lone 
valences and assuming standard pKa values at pH= 7. 
Next, conjugate gradient minimization was applied for 
3A steps to remove bad contacts introduced by proto- 
nation [N is the total number of protein atoms). A de- 
tergent micelle was then created. The micelle was made 
by placing on a circle of 8 LDAO molecules in the first 
quadrant at the z = plane, with the heads pointing to 
the exterior and the tails pointing to the origin. Sym- 
metry transformations were applied about x and y axes 
to make a ring of 32 LDAO out of the first quadrant 
LDAO molecules. The ring was then copied and trans- 
lated along the z-axis to create four more rings, each 7 A 
apart. The protein was then rotated to align the quasi- 
C2 axis with the z-axis, and translated so that the origin 
overlapped with the protein's center of mass. The rings 
were placed around the protein making a tight fit, which 
covered almost the entire alpha helix region. 

To help to form a micelle, the protein was allowed to 
relax by conjugate gradient minimization for another 3^^ 
steps while the LDAO were kept in place using a har- 
monic positional restraining force of 200 kcal/A^. Then, 
the force was removed and the system was slowly heated 
in a vacuum at a rate of 20 K/ps until 200 K. After 
this short time, the LDAO shell melted and collapsed 
into a tight micellar structure around the reaction cen- 
ter complex. The total energy of the system, the sum 
of the van dcr Waals and electrostatic energy, decreased 
during heating by several thousand kcal/mol indicating 
the creation of a more stable structure. This step was 
different from Ceccarelli and Marchi's approach^ which 
required heating to 400 K for several hundred ps in order 
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to form a compact, equilibrated micelle. Once the mi- 
celle was formed, the system charge was neutralized with 
6 sodium ions, while another 30 NaCl pairs were added 
to keep the system at an approximate 0.15 M salt con- 
centration. Then, a total of 10,506 waters (10,503 in the 
charge-separated state) were added to form a truncated 
octahedron of the simulation cell. 

Following the addition of the solvent and counterions, 
each system was run through an additional equilibration 
procedure. First, water was allowed to relax along a con- 
jugate gradient minimization for 3N steps, while the mi- 
cellar protein was held fixed with a weaker restraining 
force of 10 kcal/ A^. Next, the full solvated micellar sys- 
tem was allowed to relax for another 3A^ steps to remove 
any remaining bad contacts. Following this minimiza- 
tion, the solvated system was heated again from K to 
the desired temperature for 30 ps (NVT ensemble). Af- 
ter temperature equilibration, the volume was allowed to 
expand in a 2 ns NPT run, which stabilized in less than 
200 ps. Once the density was equilibrated, NPT produc- 
tion runs lasting 5-15 ns (1-5 ns at T = 77 K) were used 
to calculate the averages. 

A single 2 fs timestep for all MD simulations was em- 
ployed, and SHAKEi^^ was used to constrain covalent 
bonds to hydrogen atoms. For constant temperature and 
pressure, the system was coupled to a Berendsen thermo- 
stat and barostat, respectively. The long-range electro- 
static interactions were handled using a smooth particle 
mesh Ewald summation with a 10 A limit in the direct 
space sumJ^ 



APPENDIX B: ATOMIC CHARGES AND 
CHARGE TRANSFER WITHIN THE SPECIAL 
PAIR 

The partial charges of the electron transfer cofactors 
arc not provided by the Amber force field and need to 
be taken from quantum calculations. Due to the large 
size of the molecules, we modified the bacteriochlorophyll 
(Bchl) cofactors by replacing their phytyl side chains with 
methyl groups. The quantum calculations of these mod- 
ified molecules were performed using GAMESS(US)^^ 
(B3LYP DFT/3-21G) and converted to partial charges 
by CHELPG protocol, also implemented in GAMESS. 
The charge distribution of atoms of the phytyl chain was 
assumed to be the same for the neutral and charged co- 
factors and was calculated using the Antechamber mod- 
ule from Amber which employs the empirical AMl-BCC 
method. The full sets of atomic charges (with phytyl 
chains) given in Table S2 (supporting information) was 
used in the MD simulations. Similarly, the distribution 
of charge in the final charge-transfer state was obtained 
from DFT partial charges of a negatively charged Bchl~ 
anion radical and a positively charged cation radical 
Bchl+ (Table S2). The positive charge was distributed 
unequally between the two Bchls of the special pair, with 
2/3 of the positive charge residing on the L subunit (Pl) 



and 1/3 of the positive charge residing on M subunit 
(Pm), as suggested by ENDOR studies of Rhodobac- 
ter sphaeroidesr^^ The set of Aq^ charges {k runs over 
the atoms of the cofactors) obtained by subtracting the 
atomic charges in the initial neutral state from the ion- 
ized state were used to calculate the Coulomb part of the 
donor-acceptor energy gap. The partial charges on the 
protein atoms were taken from the Amber FF03 force 
field^ and the TIP3P force fieldii^ was used for the 
partial charges of water. 

We used Stark spectroscopy data by Lockhart and 
Boxer— as the starting point for determining the param- 
eters of the charge-transfer state of the photoexcited spe- 
cial pair. The change in the absorption dipole moment 
within the special pair is about /cA^ = 7 D larger than 
in an isolated bacteriochlorophyll, where /c ~ 1.2 is the 
cavity field correction factor. If this change of the dipole 
moment difference, measured at 77 K, is connected to 
the mixing between the covalent (Pa/-Pl)* and charge- 
separated, (Pm'Pl)* states of the exited special pair, 
then this change in the dipole moment can be written as 



(Bl) 



where AjUct is the dipole moment of the fully ionized 
state (Pm-P^)* and ncx is the population of this state 
at the given energy gap between the neutral and ionized 
states. In terms of the two-state Hamiltonian, this pop- 
ulation is given as 
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where Ae is the difference between diabatic energies of 
the neutral and ionized states (eq [37]) and Aujp is the 
adiabatic energy gap between the eigenvalues of the two- 
state Hamiltonian 
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Here, J is the electronic coupling clement between the 
neutral and ionized states of the excited special pair and, 
following Lathrop and Friesner^^ we assume that the 
charge-transfer state (Pj^-P^)* is predominantly mixed 
with the lower excitonic state of the dimer. 

When two bacteriochlorophyll radicals, P^ and P^^ are 
placed at their crystallographic positions, the resulting 
dipole moment of the fully ionized state is A/icT = 40.2 
D. This implies that average charge mixing between the 
two states at 77 K is ncT(77K) = 0.143. In order to 
determine the coupling parameter J from this number 
we used the model vibronic Hamiltonian of Friesner and 
co-workers which was shown to reproduce a number of 
experimental properties (absorption, circular dichroism, 
polarized absorption) In this model, the difference of 
energies between the ionized charge-transfer and neutral 
states of the special pair is 2800 cm~^, which, combined 
with the population of charge-transfer state, gives J = 
979 cm^^ and Ae = 1998 cm^^. The former value falls in 
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between 600 cm^^ used by Lathrop and Friesner— and 
1450 cm~^ used by Rengeri^ 

The electronic mixing between the neutral and ionized 
states of the special pair will make its excited state more 
polarizable than the ground state. The change in the po- 
larizability associated with charge transfer can be readily 
calculated from the two-state polarizability model which 
gives 

Aa = 2A^ll^,J^/{Aef (B4) 

With the parameters calculated above, this equation 
gives Aa = 707 A^, consistent with Aa = 460 - 745 
A'^ reported from fitting the Stark spectra.— 

The energy gap Ae was obtained in by fitting the spec- 
tra at 77 and is not directly suitable for our simula- 
tions at higher temperatures. The average energy gap be- 
tween neutral and ionized states is made by the gas-phase 
gap Ae^^^ and a shift by polar and induction interactions 
with the protein/water solvent feg I37|l . In order to ex- 
tract this shift, we have run a short (1 ns) MD simula- 
tion of the reaction center at 77 K from which the solvent 
shift was determined to be —0.974 eV. This number al- 
lowed us to determine the gas-phase gap of Ae^'^^ = 1.222 
eV which was used in the simulations of the polarizable 
special pair. The simulations required modification of 
Sander module of AMBER such that the instantaneous 
energy gap and special pair charges arc recalculated at 
each fifth time step of the MD run. 

APPENDIX C: FITTING EXPERIMENTAL 
KINETIC DATA 

Our model was applied to two sets of experimental 
kinetic data, the temperature dependence of the pri- 
mary rate from Fleming et al^ and time-resolved decays 
of the population of the photoexcited special pair from 
Wang et al^ For the latter set of data, recordered at 
T = 300 K, wc used the solvent reorganization energy 
from our MD simulations with the non-crgodic correction 
extracted from the Stokes shift dynamics fea Wf\ . Pop- 
ulation decays were calculated by solving the diffusion- 
reaction Fokker-Planck equations (egs HSHS^ for the mu- 
tants used in the experiment (Table SI). In contrast to 
Wang et al^ who used three fitting parameters in their 
analysis, the reorganization energy is fixed here by MD 
simulations and only the gas-phase energy gap AE^'^^ 
and the electron transfer matrix clement V were varied 
to fit the rate constant at 300 K for mutant L170ND, 
which is very close to the wild type reaction center with 
the difference between their mid-point potentials of only 
-0.007 eV.i^ This fit has resuffed in AE^'''' = 1.86 eV 
and V = 41.5 c m~^ (> 60 cm~^ was identified for this 
parameter in ref Il44l ). 

Special pair mutants introduce electrostatic perturba- 
tions at the location of two sandwiched bacteriochloro- 
phylls without significantly affecting the solvation com- 
ponent of the reaction Gibbs energy. This implies the 
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FIG. 17: Correlation between the change of the inho- 
mogeneous Coulomb energy gap AAi5i*^i-j between mu- 
tants and wild-type reaction center and the corresponding 
changes in the midpoint redox potential IS.IS.Em reported 
experimentally.— The values of AAi5i*J^ij are obtained from 
fitting the theoretical curves for the population decay to 
experimenl4^ (Figure 1121) while keeping the electron trans- 
fer matrix element and the gas-phase gap constant. The data 
shown in the plot are also listed in Table SI (supporting in- 
formation). The dashed line indicates the unitary slope to 
guide the eye. 



variation of the inhomogeneous Coulomb component of 
the energy gap, Ai?Pj^. This component was varied in fit- 
ting the experimental P(t^ curves of other mutants while 
keeping the gas-phase gap and the electron transfer ma- 
trix element constant. The variation of AE^^^ with muta- 
tion relative to the wild-type reaction center then closely 
follows changes in the midpoint potential (Figure fT7|) . 

The gas-phase shift AE^^^ and the matrix element V ^ 
obtained from fitting the rate at 300 K, were then used to 
calculate the temperature dependence of the rate. This 
calculation is complicated by the fact that all solvation 
energies and solvation relaxation times depend on tem- 
perature. The MD simulations do not provide sufficient 
accuracy to reliably estimate the temperature change 
of the Coulomb part of the reorganization parameters. 
Their temperature dependence was estimated from lin- 
ear interpolations with the slopes listed in Table IIVI (see 
also the discussion in sec IIII D[) . The temperature de- 
pendence of the induction component of the average gap 
is the main ingredient in reproducing the slope of the 
rate correctly. This component is converged exception- 
ally well in MD simulations which were used to produce 

A£;'"d(r). 

The temperature dependence of the Stokes shift re- 
laxation time makes the non-ergodicity correction fac- 
tor /'^'(/cet) temperature-dependent as well. It turned 
out that only the relaxation time T2{T) is strongly 
temperature-dependent in the Stokes shift correlation 
function approximated by three components according 
to eg 1211 The first two components were given constant 
values, Tfj = 0.1 ps and ti — 2.5 ps, and the longest re- 
laxation time was given the Arrhenius temperature law 
lnr2(T) = 0.936 + Er/T with Er = 1212.3 K. The results 
of these calculations, shown in Figure 1141 ^rc in a good 
agreement with the data by Fleming et alX in the range 
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TABLE IV: Parameters used to fit the charge-separation rate at 300 K and to produce the temperature dependence of the 
rate." All parameters refer to the wild-type reaction center; temperature derivatives taken at 300 K are in K~^. 

Agg=^°, eV V, cm-^ dlnAE^/dT d\n AE''"^ / dT dln\^/dT dlnA'"''/rfr 
1.86 41.5 5.2 X lO'" 4.8 x 10'"^ -1.3 x 10'^ 1.1 x 10'^ 

"Temperature dependence of the Stokes shift correlation func- 
tion was produced by using the following parameters in cq 1221 
Aa =0.172, TO = 0.1, Ai = 0.063, n = 2.5 ps, t§ = 2.55 ps, 
Er = 1212 K with T2{T) = r| cxp[l3Er]. 
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